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i In the Color Glass Condensate (CGC) effective field theory, when two 

large sheets of Colored Glass collide, as in a central nucleus-nucleus col- 
lision, they form a strongly interacting, non-equilibrium state of matter 
called the Glasma. How Colored Glass shatters to form the Glasma, the 
C^hI properties of the Glasma, and the complex dynamics transforming the 

1^ ■ Glasma to a thermalized Quark Gluon Plasma (QGP) are questions of 

central interest in understanding the properties of the strongly interacting 
matter produced in heavy ion collisions. In the first of these lectures, we 
shall discuss how these questions may be addressed in the framework of 
. , particle production in a field theory with strong time dependent external 

■ sources. Albeit such field theories are non-perturbative even for arbitrarily 

weak coupling, moments of the multiplicity distribution can in principle 
be computed systematically in powers of the coupling constant. We will 
demonstrate that the average multiplicity can be (straightforwardly) com- 
puted to leading order in the coupling and (remarkably) to next-to-leading 
order as well. The latter are obtained from solutions of small fluctua- 
tion equations of motion with retarded boundary conditions. In the second 
lecture, we relate our formalism to results from previous 2-1-1 and 3-1-1 
dimensional numerical simulations of the Glasma fields. The latter show 
clearly that the expanding Glasma is unstable; small fluctuations in the 
initial conditions grow exponentially with the square root of the proper 
time. Whether this explosive growth of small fluctuations leads to early 
thermalization in heavy ion collisions requires at present a better under- 
standing of these fluctuations on the light cone. In the third and flnal 
lecture, motivated by recent work of Bialas and Jezabek [1] , we will discuss 
how the widely observed phenomenon of limiting fragmentation is realized 
in the CGC framework. 
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Introduction 

The theme of these lectures at the 46th course of the Zakopane school is 
multi-particle production in hadronic collisions at high energies in the Color 
Glass Condensate (CGC) effective field theory. There has been tremendous 
progress in our theoretical understanding in the seven years since one of the 
authors last lectured here. At that time, the author's lectures covered the 
state of the art (in the CGC framework) in both deeply inelastic scattering 
(DIS) studies and in hadronic multi-particle production [2]. A sign of rapid 
progress in the field is that there were several talks and lectures at this 
school covering various aspects of this physics in DIS alone. Wc will restrict 
ourselves here to developments in our understanding of multi-particle pro- 
duction in hadronic collisions. Another significant development in the last 
seven years has been the large amount of data from the Relativistic Heavy 
Ion Collider (RHIC) at BNL, key features of which were nicely summarized 
in the white-papers of the experimental collaborations [3] culminating in 
the announcement of the discovery of a "perfect fluid" at RHIC. The ex- 
citing experimental observations were discussed here in the lectures of Ja- 
cak [4] . The RHIC data have had a tremendous impact on the CGC studies 
of multi-particle production. While we will discuss RHIC phenomenology, 
and indeed specific applications of theory to data, our primary focus will be 
on attempts to develop a systematic theoretical framework in QCD to com- 
pute multi-particle production in hadronic collisions. Some applications of 
the CGC approach to RHIC phenomenology were also covered at this school 
by Kharzeev as part of his lectures [5]. For recent comprehensive reviews, 
see Ref. [6,7]. 

The collider era in high energy physics has made possible investigations 
of QCD structure at a deep level in studies of multi-particle final states. 
Much attention has been focused on the nature of multi-particle production 
in jets; for a nice review, see Refs. [8,9]. The problem is however very gen- 
eral. Theoretical developments in the last couple of decades suggest that 
semi-hard particle production in high energy hadronic interactions is domi- 
nated by interactions between partons having a small fraction x of the lon- 
gitudinal momentum of the incoming nucleons. In the Regge limit of small 
X and fixed momentum transfer squared (corresponding to very large 
center of mass energies ^/s) the Balitsky-Fadin-Kuraev-Lipatov (BFKL) 
evolution equation [10,11] predicts that parton densities grow very rapidly 
with decreasing x. A rapid growth with x, in the gluon distribution, for 
fixed ^ ^QCD' observed in the HERA experiments [12,13]. (It is 
not clear however that the observed growth of the gluon distribution is a con- 
sequence of BFKL dynamics [14].) Because the rapid growth in the Regge 
limit corresponds to very large phase space densities of partons in hadronic 
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wave functions, saturation effects may play an important role in hadronic 
collisions at very high energies [15-18]. These slow down the growth of 
parton densities relative to that of BFKL evolution and may provide the 
mechanism for the unitarization of cross-sections at high energies. 

The large parton phase space density suggests that small x partons can 
be described by a classical color field rather than as particles [19-21]. Light 
cone kinematics (more simply, time dilation) further indicates that there is 
a natural separation in time scales, whereby the small x partons are the 
dynamical degrees of freedom and the large x partons act as static color 
sources for the classical field. In the McLerran-Venugopalan (MV) model, 
the large color charge density of sources is given by the density of large 
X partons in a big nucleus which contains 3^ valence quarks (where A is 
the atomic number of the nucleus). In this limit of strong color sources, 
one has to solve the non-linear classical Yang- Mills equations to obtain the 
classical field corresponding to the small x parton modes. This procedure 
properly incorporates, at tree level, the recombination interactions that are 
responsible for gluon saturation. In the MV model, the distribution of large 
X color sources is described by a Gaussian statistical distribution [19,22]. A 
more general form of this statistical distribution, for SU (iVc) gauge theories, 
valid for large A and moderate x, is given in Refs. [23,24]. 

The separation between large x and small x, albeit natural, is somewhat 
arbitrary in the MV model; the physics should in fact be independent of 
this separation of scales. This property was exploited to derive a functional 
renormalization group (RG) equation, the JIMWLK equation, describing 
the evolution of the gauge invariant source distributions to small x [25-32]. 
The JIMWLK functional RG equation is equivalent to an infinite hierarchy 
of evolution equations describing the behavior of multi-parton correlations 
at high energies first derived by Balitsky [33]. A useful (and tremendously 
simpler) large Nc and large A mean-field approximation independently de- 
rived by Kovchegov [34], is commonly known as the Balitsky-Kovchegov 
equation. The general effective field theory framework describing the non- 
trivial behavior of multi-parton correlations at high energies is often referred 
to as the Color Glass Condensate (CGC) [35-37]. 

Several lecturers at the school discussed the state of the art in our un- 
derstanding of the small x wavefunction [38] . For previous discussions at re- 
cent schools, see Refs. [39-41]. To compute particle production in the CGC 
framework, in addition to knowing the distribution of sources in the small 
X nuclear wavefunction, one must calculate the properties of multi-particle 
production for any particular configuration of sources. In this paper, we 
will assume that the former is known. All we require is that these sources 
(as the RG equations tell us) are strong sources, parametrically of the order 
of the inverse coupling constant, and are strongly time dependent. We will 
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describe a formalism to compute multi-particle production for an arbitrary 
distribution of such sources ^. 

These lectures are organized as follows. In the first lecture, we shall 
describe the formalism for computing particle production in a field theory 
coupled to strong time-dependent external classical sources. We will con- 
sider as a toy model a 0^ scalar theory, where (f) is coupled to a strong exter- 
nal source. Although the complications of QCD - such as gauge invariance 
- are very important, many of the lessons gained from this simpler scalar 
theory apply to studies of particle production in QCD. We will demonstrate 
that there is no simple power counting in the coupling constant g for the 
probability P„ to produce n particles. A simple power counting however 
exists for moments of Pn- We will discuss how one computes the average 
multiplicity and (briefly) the variance. With regard to the former, we will 
show how it can be computed to next-to-leading order in the multiplicity. 
We will also discuss what it takes to compute the generating function for 
the multiplicity distribution to leading order in the coupling. 

In lecture II, we will relate the formal considerations developed in lecture 
I, to the results of real time numerical simulations of the average multiplic- 
ity of gluons and quarks produced in heavy ion collisions. Inclusive gluon 
production, to lowest order in the loop expansion discussed in lecture I, is 
obtained by solving the classical Yang-Mills equations for two color sources 
moving at the speed of light in opposite directions [42-44]. This problem 
has been solved numerically in [45-50] for the boost-invariant case. The 
multiplicity of quark-pairs is computed from the quark propagator in the 
background field of [45-50] - it has been studied numerically in [51,52]. A 
first computation for the boost non-invariant case has also been performed 
recently [53,54]. It was shown there that rapidity dependent fluctuations 
of the classical fields lead to the non-Abelian analog of the Wcibel insta- 
bility [55], first studied in the context of electromagnetic plasmas. Such 
instabilities may be responsible for the early thermalization required by 
phenomenological studies of heavy ion collisions. We will discuss how a 
better understanding of the small quantum fluctuations discussed in lecture 
I may provide insight into early thermalization. 

In the third and final lecture, we will discuss how the formalism outlined 
in lecture I simplifies in the case of proton-nucleus collisions. (A similar 
simplification occurs in hadron-hadron collisions at forward/backward ra- 
pidities where large x's in one hadron (small color charge density) and small 
X in the other (large color charge density) are probed.) At leading order 
in the coupling, lowest order in the proton charge density, and all orders 

^ That one can separate the properties of partons in the wave function from those in 
the final state is a statement of factorization. This has not yet been proven. We will 
briefly describe work in that direction in these lectures. 
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in the nuclear color charge density, analytical results are available for both 
inclusive gluon and quark production ^. In the former case, the analytical 
formula can be written, in k±_ factorized form, as the product of uninte- 
grated distributions in both the proton and the nucleus convoluted with the 
matrix element for the interactions squared. This formula is used exten- 
sively in the literature for phenomenological applications. We will discuss 
one such application, that of limiting fragmentation. 



Lecture I: How Colored Glass shatters to form the Glasma 

As outlined in the introduction, the Glasma is formed when two sheets 
of Colored Glass collide, producing a large number of partons. A cartoon 
depicting this collision is shown in fig. 1. In the CGC framework, it is 




Fig. 1. Typical contribution to gluon production in the collision of two sheets of 
Colored Glass. The dots denote the color sources that describe the fast partons in 
the CGC framework. 

expected that observables can be expressed as^ 

(O), = j [Dp^\ [Dp2] [Pi] [P2] 0[pi,P2] , (1) 

where they are first computed as a functional of the color charge densities pi 
and p2 of the two nuclei and then averaged over all possible configurations 
of these sources, with the likelihood of a particular configuration at a given 
rapidity Y (= ln(l/x)) specified by the weight functional _y [pi] and 

^ For quark production, k±_ factorization breaks down even at leading order in pA 
collisions [117]. 

^ As we shall see in lecture II, this formula may also require an average over some 
quantum fluctuations, because of an instability in boost invariant classical solutions 
of the Yang-Mills equations. 



6 



W^Yb +Y [Ps] respectively. Here ibeam = \ ln(s/mp) is the beam rapidity in 
a hadronic collision (mp denoting the proton mass) with the center of mass 
energy -y/s. 

The evolution of the Wy [p]'s with rapidity is described by the JIMWLK 
equation [25-32]. For small x (large Y) and/or large nuclei, the rapid growth 
of parton densities corresponds to light cone source densities pi , p2 ~ 1/5 - 
in other words, the sources are strong. Thus understanding how two sheets 
of Colored Glass shatter to produce the Glasma requires that we understand 
the nature of particle production in a field theory with strong, time depen- 
dent sources. In this lecture, we will outline the tools to systematically 
compute particle production in such theories. More details can be found in 
Refs. [56,57]. 

Field theories with strong time dependent sources are different from field 

theories in the vacuum in one key respect. The "vacuum" in the former, even 
in weak coupling, is non-trivial because it can produce particles. Specifically, 
the amplitude from the vacuum state |Oin) to a populated state [aout) is 

<aout|Oin> 7^0. (2) 

Unitarity requires that the sum over all "out" states satisfies the identity 

^k«out|Oin> ' = 1 . (3) 



We therefore conclude that 

<0out|0in>|'<l. (4) 

In other words, the probability that the vacuum stays empty is strictly 
smaller than unity. Following the conventions of [63], we can write the 
vacuum-to-vacuum transition amplitude as 

<0out|0i„) = e'^[''] , (5) 

where iV[p] compactly represents the sum of the connected vacuum- vacuum 

diagrams in the presence of the external (in our case, strong, time depen- 
dent, colored) source p. Therefore, the inequality (4) means that vacuum- 
vacuum diagrams have a non-zero imaginary part, since |(0out|0in)| = 
exp(— 2Im V[p]). In stark contrast, for a field theory without external 
sources, cq. (4) would be an equality, and the vacuum-vacuum diagrams 
would be purely real, thereby only giving a pure phase for the vacuum- 
to- vacuum amplitude in eq. (5). They therefore do not contribute to the 
probabilities for producing particles. 
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Eq. (4) tells us that one has to be more careful in field theories with exter- 
nal sources. To illustrate how particle production works in such theories, we 
shall, for simplicity, consider a real scalar field with cubic self-interactions, 
coupled to an external source. (The lessons we draw carry over straightfor- 
wardly to QCD albeit their implementation is in practice significantly more 
complex.) The Lagrangian density is 

L ^ \d^<^d^^ - ^mV' - l*^' + P</> • (6) 

Note that the coupling g in this theory has dimensions of the mass; and 
that the theory is super-renormalizable in n = 4 dimensions. The source 
densities p(x) can be envisioned as the scalar analogue of the sum of two 
source terms p{x) = p\{x) -|- Piix) corresponding respectively in the CGC 
framework to the recoil-less color currents of the two hadronic projectiles. 

Let us now consider how the pcrturbative expansion for such a theory 
looks like in weak coupling. The power of a generic simply connected dia- 
gram is given simply by 

^n,+2K-l)(^^)n,^ (7) 

where are the number of external lines, the number of loops and Up 
the number of sources. For vacuum-vacuTuri graphs, = 0. As p ^ 1/g, 
the power counting for a theory with strong sources is given entirely by an 
expansion in the number of loops. In particular, at tree level (n^ = 0), the 
vacuum-vacuum graphs are all of order 1/ g^. The tree graphs contributing 
to the connected vacuum-vacuum amplitude in eq. (4) can be represented 
as 

^V\j]^^Y,V= 1^ .i]m[ ^k}^ ^ ■■ (8) 

conn 

There are also loop contributions in this expression which we have not rep- 
resented here. 

To proceed with the perturbative computation, we need to consider the 
analog of the well known Cutkosky rules for this case. For each diagram in 
the computation, begin by assigning for each vertex and source, two kinds 
of vertices denoted by -I- or — . A vertex of type -|- is the ordinary vertex and 
appears with a factor —ig in Feynman diagrams. A vertex of type — is the 
opposite^ of a -|- vertex, and its Feynman rule is +ig. Likewise, for insertions 
of the source p, insertions of type + appear with the factor +ip{x) while 



Because the coupling constant g is real in an unitary theory, the vertex of type — is 
also the complex conjugate of the vertex of type -|-. 
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insertions of type — appear instead with —ip(x). Thus for each Feynman di- 
agram iV in eq. (6), containing only + vertices and sources (denoted hence- 
forth as iVs^j^ ^}) contributing to the sum of connected vacuum- vacuum 

diagrams, we obtain a corresponding set of diagrams by assigning the 

symbol = it to the vertex i of the original diagram (and connecting a ver- 
tex of type e to a vertex of type e' with a propagator G^^, - to be discussed 
further shortly). 

The generalized set of diagrams therefore includes 2" such diagrams if 
the original diagram had n vertices and sources. Using recursively the so- 
called "largest time equation" [58,59], one obtains the identity, 

+ = -2 Im F = - ^ , (9) 

where the prime in the sum means that we sum over all the combinations 
of ej's, except the two terms where the vertices and sources are all of type 
-|- or all of type — . (There are therefore 2" — 2 terms in this sum.) 

We now need to specify the propagators connecting the it vertices and 
sources. The usual Feynman (time-ordered) free propagator is the propa- 
gator connecting two vertices of type -|-, i.e. G!|__,_. It can be decomposed 
as 

G\+{x,y) ^ e{x^ - yO) G\{x,y) + ^?(yO - x^) Gl_{x,y) , (10) 

which defines the propagators G^_^ and G^_. Likewise, the anti-time- 
ordered free propagator G° _ is defined as ^ 

G^__ix,y) ^ eix^ - yO) Gl_ix,y) + e{y^ - x^) +(x,y) . (11) 

The Fourier transforms of the free propagators G^^, for our scalar theory 

are 

Gl^{p) = - G^__{p)= , . , 

G\{p) = 27r^(pO)(5(p2 - m^) , G^_(p) = 27r^( V)5(p^ - m^) ■ (12) 

For a given term in the right hand side of eq. (9), one can divide the 
diagram in several regions, each containing only + or only — vertices and 
sources. (There is at least one external source in each of these regions be- 
cause of energy conservation constraints.) The -|- regions and — regions of 

® The notations for the propagators are those that appear in the Schwinger-Keldysh 
formahsm developed initially for field theories at finite temperature (see [64,65]). This 
identification, as we shall see later, is not accidental. The propagators so defined are 
not independent; they are related by the identity G° _|_ -|- G° _ = G° + + G° _. 
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the diagram are separated by a "cut" , and one thus obtains a "cut vacuum- 
vacuum diagram". At tree level, the first terms generated by these cut- 
ting rules (applied to compute the imaginary part of the sum of connected 
vacuum- vacuum diagrams in eq. (9)) are 



+ ■■■ 



(13) 



The + and — signs adjacent to the grey line in each diagram here indicate the 
side on which the set of -|- and — vertices is located. As one can see, there are 
cuts intercepting more than one propagator. The sum of the diagrams with 
r cut propagators is denoted by br/g^ - the identity in eq. (13) (and eq. (9)) 
is a statement of unitarity. These br are sometimes called "combinants" in 
the literature [60]. 

It is important to note that cut connected vacuum-vacuum diagrams 
would be zero in the vacuum because energy cannot flow from one side of 
the cut to the other in the absence of external sources. This is of course 
consistent with a pTire phase in eq. (9). This constraint on the energy flow 
is removed if the fields are coupled to time-dependent external sources. 
Cut vacuum- vacuum diagrams, and therefore the imaginary part of vacuum- 
vacuum diagrams, differ from zero in this case. 

We now turn to the probabilities for producing n particles. The proba- 
bility to produce one particle from the vacuum can be parameterized as 

Pi = e"i^^'-'- ^, (14) 

r 

where 6i, a series in (n > 0) is obtained by summing the 1-particle cuts 
through connected vacuum-vacuum diagrams. The exponential prefactor is 
the square of the sum of all the vacuum- vacuum diagrams, which arises in 
any transition probability. The probability P2 for producing two particles 
from the vacuum contains two pieces. One is obtained by squaring the bi/g^ 
piece of the probability for producing one particle (dividing by 2 for identical 
particles) - in this case, the two particles are produced independently from 
one another. The other 62/5^ is a "correlated" contribution from a 2-particle 
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cut through connected vacuum-vacuum diagrams. We therefore obtain 



P2 = e~^^'-' 



2! 5 



9^ 



(15) 



In a similar vein, the probability P3 can be shown to consist of three pieces. 
One ( "uncorrelated" ) term is the cube of 61/5^ (preceded by a symmetry 
factor 1/3!). Another is the combination 6162/3^, corresponding to the 
production of two particles in the same subdiagram with the third produced 
independently from the first two. Finally, there is the "correlated" three 
particle production probability 63/5^ corresponding to the production of 
three particles from the same diagram. The sum of these three pieces is 
thus 

Ibf 6162 63- 

Some of the graphs contributing to 61, 62 and 63 are shown in fig. 2. 
Following this line of inductive reasoning, one obtains a general formula 



(16) 



^= ^Hr ^-h 
+ ••• 

+ ••• 

+ ••• 

Fig. 2. Examples of cut diagrams contributing to foi, 62 and 63. 
for the production of n particles 

p=0 aiH \-ap=n " 

for any n. In this formula, p is the number of disconnected subdiagrams 
producing the n particles, and br/g^ denotes the sum of all r-particle cuts 
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through the connected vacuum-vacuum diagrams. This formula gives the 
probabiUty of producing n particles to all orders in the coupling 5 in a 
field theory with strong external sources. It is important to realize that all 
the details of the dynamics of the theory under consideration are hidden 
in the numbers 6^) and that eq. (17) is a generic form for transition proba- 
bilities when many disconnected graphs as well as vacuum-vacuum graphs 
contribute. This formula is therefore equally valid for QCD. 

It is useful to introduce a generating function for these probabilities, 

00 ^ 
F{z) = ^ P„ = exp [- ^ br (/■ - 1) ] . (18) 

n=l ^ r 

One can use this object in order to compute moments of the distribution of 
probabilities 



{n) = {lnFy{z = l) = ^J2^br. 

in') - {nf = {lnF)"{z = 1) = 1 6, , (19) 

9 J, 

where each "prime" denotes a derivative with respect to z. Note that 
F{z = 1) = J2n-^n = 1- This demonstrates explicitly that the exponential 
prefactor in eq. (17) is essential for the theory to be unitary. Though we 
derived eqs. (17) and (18) independently, we were alerted by Dremin [61] 
that an earlier version of the formulas in eqs. (17) and (18) was derived by 
Gyulassy and Kauffmann [60] nearly 30 years ago also using general combi- 
natoric arguments that did not rely on specific dynamical assumptions. 

These combinatoric rules for computing probabilities (and moments 
thereof) in field theory with strong external sources can be mapped on 
to the AGK cutting rules derived in the context of reggeon field theory [62] 
by writing eq. (17) as 

n 

Pn = Y.Pni. (20) 
p=0 

(c) 

where Pn,p denotes the probability of producing n particles in p cut sub- 
diagrams. One can ask directly what the probability of p cut sub-diagrams 
is by summing over n to obtain 

n=p y \ y / 
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This is a Poisson distribution, which is unsurprising in our framework, be- 
cause disconnected vacuum-vacuum graphs are uncorrelated by definition. 
The average number of such cut diagrams is simply 

+00 ^ 

(ncut) = ^p7^p = — ^ 6^ . (22) 

p=0 ^ r 

An exact identification with Ref. [62] is obtained by expanding the expo- 
nential in eq. (21) to order m ~ p, and defining 



im-p)\p\ [ ) V 5^ ' ' ^ ' 



where TZp^m is the probability of having p cut sub-diagrams out of m sub- 
diagrams (with m — p being the number of uncut diagrams). This distribu- 
tion of probabilities can be checked to satisfy the relations 



if m > 2 , ^^•7^p,^ = 0, 



p=i 



ifm>3, y^p(p - l)^p,m = , 

p=2 



(24) 



This set of identities is strictly equivalent to the eqs. (24) of Ref. [62] where 
m — p and p there are identified as the numbers of uncut and cut reggeons 
respectively. The first relation means that diagrams with two or more sub- 
diagrams cancel in the calculation of the multiplicity. These relations are 
therefore a straightforward consequence of the fact that the distribution of 
the numbers of cut subdiagrams is a Poisson distribution. They do not 
depend at all on whether these subdiagrams are "reggeons" or not. In the 
AGK approach, the first identity in eq. (24) suggests that the average num- 
ber of cut reggeons (ricut) can be computed from diagrams with one cut 
reggeon alone. The average multiplicity satisfies the relation 

(n) = (ncut) {n)i , (25) 

where (n)i is the average number of particles in one cut reggeon. Computing 
this last quantity of course requires a microscopic model of what a reggeon 
is. 

Before going on, it is useful to summarize what we have learnt at this 
stage about field theories with strong external sources. We derived a general 
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formula in eq. (17) for the probability to produce n particles in terms of cut 
connected vacuum- vacuum diagrams, where br is the sum of the terms with 
r cuts. This formula is a purely combinatoric expression; it does not rely on 
the microscopic dynamics generating the br- Nevertheless, it tells us several 
things that were not obvious. Firstly, the probability distribution in eq. (17) 
is not a Poisson distribution, even at tree level, if any 6^ 7^ for r > 1. It is 
often assumed that classical dynamics is Poissonian. We see here that the 
non-trivial correlations (symbolized by non-zero br terms with r > 1 ) in 
theories with self-interacting fields can produce significant modifications of 
the Poisson distribution. 

Another immediate observation is that even the probability to produce 
one particle (eq. (14)) is completely non-perturbative in the coupling con- 
stant g for arbitrarily small coupling. In other words. Pi cannot be ex- 
pressed as an analytic expansion in powers of g. Therefore, while weak cou- 
pling techniques are certainly valid, such theories (the CGC for instance) 
are always non-perturbative. Interestingly, we will see shortly that a simple 
expansion in powers of the coupling exists for moments of the probability 
distribution. Finally, we saw that there was a simple mapping between 
the cutting rules first discussed in Ref. [62] and those for cut connected 
vacuum-vacuum graphs in field theories with strong sources. 

In the rest of this lecture, we shall sketch the derivation of explicit ex- 
pressions for the n-particle probabilities and for the first moment of the 
multiplicity distribution. Specifically, we will outline an algorithm to com- 
pute the average multiplicity up to next-to-leading order in the coupling 
constant. The probability for producing n particles is given by the expres- 
sion 

KPi---p„out|Oi„>|' , (26) 

where Ei = \/ Pi^ + m^. The well known Lehman-Symanzik-Zimmerman 
(LSZ) reduction formula [63] relates the transition amplitude for producing 
n particles from the vacuum to the residue of the multiple poles of Green's 
functions of the interacting theory. It can be expressed as 



n 



l\ (27r)32^, 



Zn/2 



i5p{xi 



iV[p] 



(27) 

where the factors of Z correspond to self-energy corrections^ on the cut 
propagators of the vacuum-vacuum diagrams. Substituting the r.h.s. of 



More precisely, they are the wavefunction renormahzation factors. 



14 



this equation into eq. (26), and noting that 

is the Fourier transform of the propagator given in eq. (12), we can write 
the probabihty P„ directly as 

P„ = ^P'^[p+,p_] e^^[^+l e-^^*!''-] , (29) 

n! p+=p-=p 

where is the operator 

V = Jdxdy ZG^_{x,y) ^J^^ " ^^O) 

The sources in the amphtude and the complex conjugate amplitude are 
labeled as /?+ and p_ respectively to ensure that the functional derivatives 

act only on one of the two factors. 
A useful and interesting identity is 

where iVgj^[p+, p-] is the sum of all connected vacuum- vacuum diagrams in 
the Schwinger-Keldysh formalism [64,65], with the source p+ on the upper 
branch of the contour and likewise, p- on the lower branch. When p+ = 
P- = p, it is well known that this sum of all such connected vacuum-vacuum 
diagrams is zero. The generating function F{z), from eqs. (18) and (29) is 
simply 

F{z) = e'^ e^^^P+'^ e-'^*^P-^ . (32) 

p+=p-=p 

Prom the expression for the operator D in eq. (30), it is clear that F{z) can 
be formally obtained by substituting the off-diagonal propagators _|_ — 
z G^^^ in the usual cut vacuum- vacuum diagrams. 

We shall now proceed to discuss how one computes the average mul- 
tiplicity {F'{z = 1)) of produced particles. Prom eqs. (32) and (30), we 
obtain 

<n) = / d^xd^y ZG'_,_{x,y) [rW(x)r(-)(y) + r(+-)(x, y)! 

(33) 

where r*^^^ and T^'^ ^ arc the 1- and 2-point amputated Green's functions 
in the Schwinger-Keldysh formalism: 

r(±)(ar) = + ^i^sAP+^P-] 
Z 5p±{x) 

r(+-VT °^+"^' By + m^ d^iV,^[p+,P-] 

^ Z Z Sp^{x)6p^{y) ■ ^'"^ 
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Diagrammatically, (n) can be represented as 

Unlike the probabilities, there is a well defined power counting for the mo- 
ments of the multiplicity distribution. This is simply because the overall 
"absorption factor" exp(— hr/g^) present in the computation each prob- 
ability, cancels when one computes averaged quantities. This is a crucial 
simplification, because it means that the moments of the distribution have 
a sensible^ perturbative expansion as a series in powers of . 

At leading order in the coupling constant, O {g~^ {g p)"-) , only the left 
diagram in eq. (35) contributes. The right diagram, that contains the con- 
nected 2-point function, is a one loop diagram; in our power counting (see 
eq. (7)) it starts at order 0{g^{gp)^). The lowest order, where we need only 
tree level diagrams, can therefore be expressed as 

R„=i: V4{ . (36) 

where the sum is over all the tree diagrams on the left and on the right of the 
propagator (represented in boldface) as well as a sum over the labels 
+/— of the vertices whose type is not written explicitly. At this order, the 
mass in G^_ is simply the bare mass, and Z = 1. 

The diagrams in eq. (36) can be computed using the Cutkosky rules 
we discussed previously. Beginning with one of the "leaves" of the tree 
(attached to the rest of the diagram by a + vertex for instance), one has 
two contributions and —G^_ for the propagators connecting it to 

the vertex just below. (The source can be factored out, because we set 
p+ = p- = p.) This difference in the propagators gives 

G° + - G° _ = Gl , (37) 

where is the free retarded propagator. (Likewise, G^^ — G^__ = G^.) 
Repeating this procedure recursively, propagators from all the "leaves" down 
to the root are converted into retarded propagators. It is well known that 
the retarded solution (f)c{x) of the classical equations of motion with the 
initial conditions lim.j.o^_^(pc{x) = and lini^o^_^d^(pc{x) = can be 
expressed as a sum of tree diagrams built with retarded propagators. The 



The usual caveats, about the convergence of such series and issues related to their 
Borel summability, apply here as well. 
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sum over all the trees on each side of the cut in eq. (36) can therefore be 
identified as 



+/- 





(38) 



Prom this discussion and eq. (33) , the leading order inclusive multiplicity 
can be expressed as 



(39) 



Using the identity e*^''^(9o + Ep)(f)c{x) = Oq (e*^'^[9o — iEp]4>c{x)) and the 
boundary conditions obeyed by the retarded classical field <j)c{x), one obtains 



E, 
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lim / d^x e^f^ [d^o - iEp] ^c{x) 



(40) 



The corresponding formula for gluon production in heavy ion collisions in 
the Color Glass Condensate framework is 



E: 



d{n) 



LO 



^ lim / d^x d^y e^^"^^-^') {d^o - iEp) {dyO + iEp) 



^ d^p IGtt^ 

phys. A 



(41) 



where is the polarization vector for the produced gluon. This is precisely 
the expression that was computed in previous real time numerical simula- 
tions of Yang-Mills equations for each configuration of color sources in each 
of the nuclei. To compute the distribution of gluons, we need to average 
over the distribution over all possible color sources as specified in eq. (1). 
We will discuss results from these simulations further in Lecture II. 

The leading order result in cqs. (40) and (41) is well known. Wc shall 
now discuss the computation to next-to-leading order in the coupling - 
to order 0{g'^{gj)^). At this order, both terms in eq. (33) contribute to 
the multiplicity. The right diagram in eq. (35) contributes with the blob 
evaluated at tree level, 



(42) 
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This contribution to the inclusive multiphcity is analogous to that of quark- 
anti-quark pair production or gluon pair production to the respective av- 
erage multiplicities for these quantities. The left diagram in cq. (35), at 
this order, contains 1-loop corrections to diagrams of the kind displayed in 
eq. (36) . A blob on one side of the cut in eq. (35) is evaluated at the 1 loop 
level (corresponding to the contribution from one loop correction to the 
classical field) while the other blob is evaluated at tree level (corresponding 
to the contribution from the classical field itself). This can be represented 
as 



(43) 



The inclusive multiplicity at NLO includes contributions from both eqs. (42) 
and (43). 

To evaluate the diagram in eq. (42), one needs to compute the propagator 

G_| {x,y) in the presence of the background field (pc- This can be done by 

solving a Lippmann-Schwinger equation for G+_ [51,56,66]. In practice, 
numerical solutions of this equation can be obtained only for retarded or 
advanced Green's functions in the background field. It turns out that one 
can express in terms of these as 

G+. = G^Gl-'Gl_G'^-'G,, (44) 

where 

Gu - Gl + GlT.Gl , 

^ G° + GO r,GO . (45) 

Here G^ (G^) is the free retarded (advanced) propagator and (T^) 
is the retarded (advanced) scattering T-matrix. Substituting eq. (45) in 
eq. (44) and using the resulting expression in the second term of eq. (33), 
the contribution of this term to the NLO multiplicity can be expressed as 



(46) 



^'^/nlo J {2Trf2Ej (27r)32E, ' ' 
One can then show that [56] 

T« {p, -q) = lim / d^x e^f [d,, - iEp] 7^g{x) , (47) 

where r]q{x) is a small fluctuation field about (j)c{x) and is the retarded 
solution of the partial differential equation 

{n + m^ + gMx)hq{x) = , (48) 
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with the initial condition r]q{x) = e^^'^ when xq — oo. Note here that 
g has the dimension of a mass. Note also that, despite being similar, the 
equation for 77 is not the classical equation of motion but is instead the 
equation of motion of a small fluctuation. This NLO contribution to the 
inclusive multiplicity can be computed by solving an initial value problem 
with boundary conditions set at x° ^ — cxd. 

The other contribution of order 0{g^{gjy'^) to the average multiplicity 
is from the diagram in eq. (43). This contribution can be written as 



(2) 
NLO 



lim f d^x e'P-'' \do - lEp] 4>c 1 (a;) 

a;o-++oo J ■' ' . 



+ C.C. 



(49) 



The one loop contribution to the classical field 



</'c,i(a;) 



E 

+/- 



X 




5 


+ 

t 





o 



(50) 



includes arbitrary insertions of the background field (t)c{x)- Following the 
discussion before eq. (37) of the Cutkosky rules in this case, it can be written 
as [56] 

4^c,i{x) = -ig J d'^y Gj^{x,y) G++{y,y) . (51) 

We have used here the identity G++{x,x) = G (a;, a;). In practice, (t>c,i{x) 

can also be obtained as the retarded solution to the equation 



{U + rr? + g(f)c{x))(f)c,i{x) = -gG++{x,x) , 



(52) 



with an initial condition such that 0c, 1 and its derivatives vanish at xq = 
— CO. The source term in this equation can be rewritten as G++(x,a;) = 
\Gr{G\)-^G%{G\)-^Ga, where G%{p) = 27r5(p2 - m^). After a little 
algebra [56], one can show that 

G++(x,x) = IJ^ ^^S{q' - m') r?(+)(x)r?(-)(x) . (53) 
Here rjq^\x) and r)g~\x) are solutions of eq. (48) with plane wave initial 



conditions at xq — —00 of r]q^\x) = e**'^ and rjg '{x) = e '^'^ respectively. 
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We observe that G++(a;,x) contains ultraviolet divergences that arise from 
the integration over the momentum q in eq. (53). They can be identified 
with the usual 1-loop ultraviolet divergences of the (p"^ field theory in the 
vacuum and must be subtracted systematically in order to obtain a finite 
result. 

To summarize, the two NLO contributions to the inclusive multiplic- 
ity, eqs. (46) and (49) can be computed systematically as follows. One 
first computes the lowest order classical field (pc{x) by solving the classical 
equations of motion, as a function of time, with the retarded boundary con- 
dition (l)c{x) = at = — oo. This computation was performed previously 
in the CGC framework [45-50]. The small fiuctuation equation of motion 
in eq. (48) is then solved in the background of (f)c{x), also with retarded 
boundary conditions at x^ = — oo for the small fluctuation field 'r]q{x). This 
is then sufficient, from eqs. (47) and (46), to compute one contribution to 
the NLO multiplicity. To compute the other, solutions of the small fiuctua- 
tion equations of motion can also be used, following eq. (53), to determine 
G++(.x,a;). Subsequent to this determination, the temporal evolution of 
the one loop classical field can be computed by solving eq. (52), again with 
retarded boundary conditions at x^ = — oo. Finally, this result can be sub- 
stituted in eq. (49) in order to compute the second contribution to the NLO 
multiplicity. 

Albeit involved and technically challenging, the algorithm we have out- 
lined is straightforward. The extension to the QCD case can be done. In- 
deed, this computation is similar to a numerical computation (performed 
by Gelis, Kajantie and Lappi [51]) of the number of produced quark pairs 
in the classical background field of two nuclei. 

An interesting question wc shall briefly consider now is whether we can 
directly compute the generating function itself to some order in the cou- 
pling; even a leading order computation would contain a large amount of 
information. Prom eqs. (32) and (30), we obtain^ 

^ = J d^xA^G°_(x,2/)[r(+)(^|x)r(-)(z|y) + r(+-)(z|a;,2/)] , (54) 

where rW(z|x) and F^'' )(z|x,y) are defined as in eq. (34), but must be 
evaluated with the substitution _|_ z G^ j. of the off-diagonal propa- 
gators. Unsurprisingly, this equation involves the same topologies as that 
for the average multiplicity in eq. (35). If we can compute the expression in 
eq. (54) even to leading order, the generating function can be determined 
directly by integration over z, since we know that F{1) = 1. 



From this relation, wc sec that the logarithm of F(z) has a well defined perturbativc 
expansion in powers of (that starts at the order g^^), while this is not the case for 
F{z) itself. 
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At leading order, as we have seen, only the first term in eq. (54) con- 
tributes and (using the same trick as in eq. (40)) eq. (54) can be written 
as 



F'(z) 



F{z) 



{2'Kf2Ep 



£x e^P-^ {d^ - lEp) ^+{z\x) 
d'y e-'P-y (5° + ^Ep)$_(^|y) 



xo=+oo 

xo=—oo 

yo=+oo 
J yo=-oo 



where ^±{z\x) correspond to the tree diagrams 



+/- 




+/- 



(55) 



(56) 



evaluated with Cutkosky's rules where the ofF-diagonal propagators 
are multiplied by a factor z. 

We will now see why computing the generating function at leading order 
is significantly more complicated than computing the average multiplicity. 
The fields in eq. (56) can equivalently be expressed as the integral equation 



\G%ix,y)[j{y)-^-^l{z\y) 



^+{z\x) = i J d^y |( 

-zG\_{x,y)[j{y)-\^l{z\y)\] 
$_(z|x) = ijd''y {zG'L+{x,y) [j{y) - ^-^l{z\y) 

-G^__{x,y)[j{y)-Ul{z\y)]} . (57) 

Note now that when = 1, $± 4>c (defined in eq. (38)) and the prop- 
agators in eqs. (57), from cq. (37), can be rearranged to involve only the 
retarded free propagator G^. It is precisely for this reason that the compu- 
tation of the inclusive multiplicity simplifies; the field (pc can be determined 
by solving an initial value problem with boundary conditions at x° — — oo. 
This simplification clearly does not occur when z ^ 1. 

In order to understand the boundary conditions for $-|-(2;|a;) in eqs. (56) 
and (57), we begin by expressing them as a sum of plane waves. 



+ 



{z\x)^ I ^J'P^ {/W(z|x°,p)e-^^-- + fi-\z\x\p)e^P-^} , 



^-i^lx)^ I {/l+)(z|x°,p)e--^-- + f[-\z\x\p)e^P-^} (58) 
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Note here that po = a/p^ + tt? is positive. Since ^±{z\x) does not obey 
the free Klein-Gordon equation, the coefficients functions must themselves 
depend on time. However, assuming that both the source p{x) and the 
coupling constant g are switched off adiabatically at large negative and 
positive times, the coefficient functions f^\z\x^,p) become constants in 
the limit of infinite time {x^ — ^ ±00). 

The technique we use for determining the boundary conditions for the 
coefficients f!^\z\x^ ,p) is reminiscent of tlic derivation of Green's theorem 
in electrostatics. We will not go into the derivation here (see Ref. [57] for 
the detailed derivation). The boundary conditions at x^ = ±00 are 

/f^(^k° = -oo,p) = 0, 
/H(z|xO = -oo,p) = 0, 
f^_^\z\x^ = +00, p) = z/|+^(z|x° = +00, p) , 
/|"^(z|x° = +00, p) = zf[~\z\x^ = +00, p) . (59) 
Using eqs. (58) and (59), we can write eq. (55) as 



F'{z) 



F(z) 



/i+)(z| + 00, p) /i-)(z| + oo,p) . (60) 



{2TTf2E. 



Therefore evaluating the generating function at leading order requires that 
we know the coefficient functions at x^ = +00. 

Unlike the case of partial differential equations with retarded boundary 
conditions, there are no straightforward algorithms for finding the solution 
with the boundary conditions listed in eq. (59). Methods for solving these 
sorts of problems are known as "relaxation processes" . A fictitious "relax- 
ation time" variable ^ is introduced and the simulation is begun at ^ = 
with functions that satisfy all the boundary conditions but not the equa- 
tion of motion. These fields evolve in ^ with the equation (preserving the 
boundary conditions for each ^) 

5^$± = (D^ + m2)$± + |$| - j{x) , (61) 

which admits solutions of the EOM as fixed points. The r.h.s. can in 
principle be replaced by any function that vanishes when $± is a solution 
of the classical EOM. This function should be chosen to ensure that the 
fixed point is attractive. A similar algorithm has been developed recently 
to study the real time non-equilibrium properties of quantum fields [67]. 

Higher moments of the multiplicity distribution can also be computed 
following the techniques described here. Interestingly, the variance (at lead- 
ing order) can be computed once one obtains the solutions of the small fluc- 
tuation equations of motion. The computation is outlined in Ref. [57] . Thus 
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both the leading order variance and the NLO inclusive multiplicity can be 
determined simultaneously. The variance contains useful information that 
can convey information about the earliest stages of a heavy ion collision. 
In particular, correlations between particles in a range of rapidity windows 
can provide insight into the early stages of a heavy ion collision [68]. This 
provides a segue for the topic of the second lecture on the properties of the 
Glasma. 

Lecture II: What does the Glasma look like and how does it 
thermalize to form a Quark Gluon Plasma ? 

In the previous lecture, we outlined a formalism to compute particle 
production in field theories with strong time dependent sources. As argued 
previously, the Color Glass Condensate is an example of such a field theory. 
In the CGC framework, the high energy factorization suggested by eq. (1) 
is assumed to compute final states. In this lecture, we will address the 
question of how one computes in practice the initial Glasma fields after 
a heavy ion collision, what the properties of these fields are and outline 
theoretical approaches to understanding how these fields may thermalize to 
form a Quark Gluon Plasma. A cartoon depicting the various stages of the 
spacetime evolution of matter in a heavy ion collision is shown in fig. 3. 



Fig. 3. Space-time development of a nucleus-nucleus collision. A goal of the Color 
Glass Condensate approach is to describe the first stage - dominated by strong 
fields - and to match it to the subsequent descriptions by kinetic theory or hydro- 
dynamics. 



t 




strong fields — > classical EOMs 
-> Z (beam axis) 



hydrodynamics 



kinetic theory 



In the CGC effective field theory, hard (large x) parton modes in each 
of the nuclei are Lorentz contracted, static sources of color charge for the 
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soft (small x) wee parton, Weizsacker-Williams modes in the nuclei. Here 
X is the longitudinal momentum fraction of partons in the colliding nuclei. 
Wee modes with x <C A~^/^ and k±_ > Aq^j-, are coherent across the longi- 
tudinal extent of the nucleus and therefore couple to a large density of color 
sources. With increasing energy, the scale separating soft and hard modes 
shifts towards smaller values of x; how this happens can be quantified by a 
Wilsonian RG [37]. In a heavy ion collision, the color current corresponding 
to the large x modes can be expressed as 

jM,a = S'^+p-{x^)S{x-) + S''-p-2{x^)S{x+) , (62) 

where the color charge densities of the two nuclei are independent 

sources of color charge on the light cone. Let us recall that x^ = {t±z)/^/2. 
The 6 functions represent the fact that Lorentz contraction has squeezed the 
nuclei to infinitesimally thin sheets. The absence of a longitudinal size scale 
ensures that the gauge fields generated by these currents will be boost- 
invariant - they are independent of the space time rapidity r] = atanh(z/t). 

The gauge fields before the collision are obtained by solving the Yang- 
Mills equations 

D^F^''' = r , 

where = 0^ + ig[Af,,.] and F^,^ = d^A^ - di^A^ + ig[A^,Au\ are the 
gauge covariant derivative and field strength tensor, respectively, in the 
fundamental representation and [A^, .] denotes a commutator. 

Before the nuclei collide {x^ < 0) , a solution of the equations of motion 
is [19,20] 

A± = , A' = e,{x-)e,{-x-^)c^-y{xA_) + ee{x+)ee{-x-)a\{xi_) , (63) 

where, here and in the following, the transverse coordinates x, y are labeled 
by the Latin index i = 1, 2. The subscript e on the ^-functions denote that 
they arc smeared by an amount e in the respective x^ light cone directions. 
We require that the functions a\y^{x±) (m = 1,2 denote the labels of the 
colliding nuclei) are such that = - they are pure gauge solutions of 
the equations of motion. The gauge fields, just as the Weizsacker-Williams 
fields in QED, are therefore plane polarized sheets of radiation before the 
collision. The functions satisfy 

-AaL = Pm{xi) . (64) 
This equation has an analytical solution given by [22,25] 

< = -- e'^- a' e-*^- , VlK^ = -gpm. (65) 
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To obtain this result one has to assume path ordering in x respectively for 
nucleus 1 and 2; we assume that the limit e — > is taken at the end of the 
calculation. 

We now introduce the proper time r = yjt"^ — z'^ = \/2x+x~ - the initial 
conditions for the evolution of the gauge field in the collision are formulated 
on the proper time surface r = 0. They are obtained [42,43] by generalizing 
the previous ansatz for the gauge field to 

A\X-,X-^,X^) = ee{x-)0e{-X+) a\{x±) + 0e{-X-)ee{x+) ai{x±) 

+6e{x~)9e{x'^) al{x~ ,x'^ ,x±) 
= ee{x-)9e{x+) P{x-,x+, x±) , (66) 

where we adopt the Fock-Schwingcr gauge condition = x'^A~ +x~A'^ = 
0. This gauge is an interpolation between the two light cone gauges = 
on the = surfaces respectively. 

The gauge fields as,/? in the forward light cone can be determined from 
the known gauge fields ai,2 of the respective nuclei before the collision by 
invoking a physical "matching condition" which requires that the Yang-Mills 
equations Dfj,F'^^ = be regular at r = 0. The (^-functions of the current 
in the Yang-Mills equations therefore have to be compensated by identical 
terms in spatial derivatives of the field strengths. Interestingly, it leads to 
the unique solution [43] 

al{x'^,x~,x_i_) = a\{x±) +al{x±) 

f3{x^,x~,x±) = [ai,i{x±),al{x^)] . (67) 

Further, the only condition on the derivatives of the fields that would lead 
to regular solutions are drP\T=o , dra'^^\T=o = 0. 

For the purpose of solving the Yang-Mills equations for a heavy-ion 
collision on a lattice, we shall work with the r, rj co-ordinates and re-express 
the initial conditions for the fields and their derivatives in terms of the fields 
and their conjugate momenta in these co-ordinates. Our gauge condition is 
A'^ = 0, and the initial conditions in eq. (67) for the functions a\ and /? at 
r = can be expressed in terms of the fields 

A{r,x±) = Ai{T,ri,x±) , 

$(r,£c_L) = x'^A'{T,r],x±)-x~A+{T,ri,x±) , (68) 

where we have made manifest the fact that these fields are boost-invariant 
- i.e. independent of r). This is a direct consequence of the assumption in 
eq. (62) that the currents are (5-function sources on the light cone. The light 
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cone Hamiltonian in = gauge, in this case of boost invariant fields, can 
be written as [46] 



K = Tr 



(69) 



Here the conjugate momenta to the fields are the chromo-electric fields 

Ei = TdrAi , En = \^. (70) 

Note that the contribution of the hard valence current does not appear ex- 
plicitly in the A"^ = Hamiltonian expressed in (r, t]) co-ordinates. The 
dependence on the color source densities is entirely contained in the de- 
pendence of the initial conditions on the source densities. Boost invariance 
simplifies the problem tremendously because the QCD Hamiltonian in this 
case is "dimensionally reduced" to a 2 -|- 1-d (QCD -|- adjoint scalar field) 
Hamiltonian. 

In terms of these Glasma fields and their conjugate momenta, the initial 
conditions in eq. (67) at r = can be rewritten as 

A'{0,x^) = a\{x±) + ai{x±) , 
EiiO,x^) = 0, 

Er^{0,x_L) =ig[a\{x^),ai{x_i)] . (71) 

The magnetic fields being defined as = ^knvF^^ , these initial condi- 
tions suggest that 5^ 7^ and Bi = 0. Note that the latter condition 
follows from the constraint on the derivatives of the gauge field that ensure 
regular solutions at r = 0. Thus one obtains the interesting results that 
the initial Glasma fields correspond to large initial longitudinal electric and 
magnetic fields {E^j, B^ ^ 0) and zero transverse electric and magnetic fields 
{Ei,Bi = 0). This is in sharp contrast to the electric and magnetic fields 
of the nuclei before the collision (the Weizsacker-Williams fields) which are 
purely transverse! Their importance was emphasized recently by Lappi and 
McLerran [69] who also coined the term "Glasma" to describe the properties 
of these fields prior to equilibration. 

An immediate consequence of these initial conditions, as noted by Khar- 
zeev, Krasnitz and Venugopalan [70], is t hat non-zero Chern-Simons charge 
can be generated in these collisions. The dynamics of the Chern-Simons 
number in nuclear collisions however differs from the standard discussion in 
two ways. Firstly, the time translational invariance of the fields is broken 
by the singularity corresponding to the collision. Secondly, due to the boost 
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invariance of the solutions, there can be no non-trivial boost invariant gauge 
transformations. This can be seen as follows. In Ref. [70], it was shown that 
the Chern-Simons charge per unit rapidity could be expressed as 



u = 



167r2 



j (fx± ^"B^ . (72) 



Because this density is manifestly invariant under rapidity dependent trans- 
formations, such transformations (which correspond to sphaleron transi- 
tions) cannot change the Chern-Simons charge. Thus sphaleron transitions 
are disallowed for boost-invariant field configurations. Eqs. (71) tell us that 
z/(r = 0) = 0; therefore the Chern-Simons charge generated in a given win- 
dow in rapidity at a time r is simply, by definition, i^(T)(77inax — ??min)- Since 
?7's of either sign are equally likely, the ensemble average {v{t)) is zero. 
However, {u{t)'^) is non zero. Its value was computed in Ref. [70]. The 
topological charge squared per unit rapidity generated for RHIC and LHC 
collisions is about 1-2 units. In contrast, estimates of the same quantity 
in a thermal plasma are one to two orders of magnitude larger. If boost 
invariance is violated (as we shall soon discuss), sphaleron transitions can 
go, and can potentially be large. This possibility, in a different formulation, 
was discussed previously by Shuryak and collaborators [71]. 

We shall now discuss the particle distributions that correspond to the 
gauge fields and their conjugate momenta in the forward light cone. Prom 
the Hamilton equations 

the Yang-Mills equations are 

r 

drAr, = TEr, , 

drEi = TDjFji + T-^Dr^Fr^i , 

drE,, = T-^DjFjr, . (74) 

They also satisfy the Gauss law constraint 

DiEi + Dr^Er, = . (75) 



These equations are non-linear and have to be solved numerically. A lattice 
discretization is convenient because it preserves gauge invariance explicitly. 
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One can write down the analogue of the well known Kogut-Susskind Hamil- 
tonian in this case and solve eq. (74) numerically on a discretized spatial^ 
lattice with the initial conditions in cq. (71). We shall not describe the 
numerical procedure here but instead refer the reader to Refs. [46,49]. 

Solving Hamilton's equations, the average gluon multiplicity can be com- 
puted using precisely the formula we discussed previously in eq. (41). The 
result in eq. (41) is the average multiplicity for a configuration of color charge 
densities in each of the nuclei. It is an average in the sense of being the first 
moment of the multiplicity distribution. This multiplicity has to be further 
averaged over the distribution of sources -v\pi\ and W„ .v\p2\, 

beam beam"*" 

as specified in cq. (1). These weight functionals have to be specified at an 
initial scale Yq in rapidity , and are then evolved to higher rapidities by 
the JIMWLK renormalization group equation [25-32] . For the purposes of 
computing the average multiplicity in central Au-Au collisions at RHIC, i.e. 
for rapidities where evolution effects a la JIMWLK are not yet important, 
the weight functionals Wy [p] are Gaussian distributions specified in the MV 
model (discussed briefly in the introduction to these lectures): 

W,^ M = exp J ' ^^^^ 

Here = g'^fj^, where g^/Jp is the color charge squared of the sources per 
unit area. The nuclei, for simplicity, are assumed to be identical nuclei. A^ is 
the only dimensionful scale (besides the nuclear radius R) in the problem. It 
is simply related, in leading order, to the nuclear saturation scale Qs by the 
expression A^ = 2TrQ'^jNcln{Qs/2L), where L is an infrared scale of order 
■^QCD- The nuclear saturation scale, performing a simple extrapolation of 
the HERA data on the gluon distribution of the proton to Au nuclei, is of 
the order Qg ~ 1.2 — 1.4 GeV at RHIC energies in several estimates [72,73]. 
For L = 0.2 GeV, this corresponds to a value A^ = 1.66 — 1.8 GeV. Clearly, 
there are logarithmic uncertainties in this estimate at least of order 10%. 
For the rest of this lecture, we will assume the Gaussian form in eq. (76) 
for the averaging over sources; modifications to account for the (very likely) 
significant effects of small x quantum evolution will have to be considered 
at LHC energies. 

In order to compute gluon number distributions, we impose the trans- 
verse Coulomb gauge • A± = to fix the gauge freedom completely. The 
result for the number distributions, averaged over the sources in eq. (76), is 



^ The proper time r is treated as a continuous variable, that can have increments as 
smaU as required to reach the desired accuracy in the solution of the equations of 
motion. 
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computed at a time r ~ 3/As to be 



1 



dN 



1 



where fn{k±/^s) is a function of the form 



^ /„(fc±/As) , 



(77) 



ai 



exp yykj_ +m?/Tcfi 
, a2 Af ln(47rA;_L/As) 



1 



(fcx/A, < 1.5) 
{k^/Ks > 1.5) 



(78) 



with ai = 0.137, m = 0.0358 A^, T^s = 0.465 A^, and 02 = 0.0087. These 
results are plotted in fig. 4 ~ they are compared to those computed in- 
dependently by Lappi [50]. The different lines in the figure correspond to 
different lattice discretizations; the differences at large kj_ therefore indicate 
the onset of lattice artifacts, which can be eliminated by going closer to the 
continuum limit (larger lattices). 




Fig. 4. Comparison of gluon transverse momentum distributions per unit area as 
a function of fc^/A^. KNV I (circles): the number defined with taken to mean 
the lattice wave number along one of the principal directions. KNV II (squares) 
and Lappi (solid line): the number defined by averaging over the entire Brillouin 
zone and with taken to mean the frequency u!{k±). 

Prom eq. (78), the number distribution at large k± has the power law 
dependence one expects in perturbative QCD at leading order. For small 
k±, the result is best fit by a massive 2-d Bose-Einstein distribution even 
though one is solving classical equations of motion! There is an interesting 
discussion in the statistical mechanics literature that suggests that such a 
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distribution may be generic for classical "glassy" systems far from equilib- 
rium [74]. Another interesting observation is that the non-perturbative real 
time dynamics of the gauge fields generates a mass scale ra which makes the 
number distributions infrared safe for finite times. Such a "plasmon mass" 
can be extracted from the single particle dispersion relation; it behaves dy- 
namically as a function of time precisely as a screening mass does [47,53]. 
This can be seen in fig. 5. As we shall discuss shortly, this plasmon mass 
can be related to the growth rate of instabilities in the Glasma. 




Fig. 5. Time evolution of the "plasmon frequency" Wpi, for fixed g^^L = 22.6 and 
lattice spacings = 0.707,0.354,0.177 (iV_L = 32,64, 128), respectively. 

The total transverse energy and number can be obtained independently, 
from the Hamiltonian density and from a gauge invariant relaxation (cool- 
ing) technique respectively. These agree with those obtained by integrating 
eq. (77) over and can be expressed as 

7ri?2 dr) ' ' ttR^ dr] ' ' ^ ' 

where = 0.27 — 0.25 and = 0.315 — 0.3 for the wide range AsR = 
50 — 167 respectively. For larger values of AgR, the functions and /^^ have 
a weak logarithmic dependence on AsR. If we assume parton-hadron duality 
and directly compare the number of gluons from eq. (79) to the number of 
hadrons measured at = in ^/s = 200 GeV/nucleon Au-Au collisions 
at RHIC, one obtains a good agreement for A^ « 2 GeV. This value is a 
little larger than the values we extracted from extrapolations of the HERA 
data; one should however keep in mind that additional contributions to 
the multiplicity of hadrons will accrue from quark and gluon production at 
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next-to-leading order [51]. If we include these contributions, as we hope to 
eventually, will be lower than this value. 

The "formation time" tj, defined as the time when the energy density 
e behaves as 1/r, is defined as tj = l/jAg, where 7 = 0.3 in the range of 
interest. The initial energy density for times t > Tf {tj ~ 0.3 fm for = 2 
GeV)is then 

1 dE^ _ 0.26 A3 

This energy density, again for A^ = 2 GeV (and g = 2), is e « 40 GeV/fm^ 
at r = Tj^. Because the energy density is ultraviolet sensitive, this number 
is probably an overestimate because the spectrum at large k± > Ag in 
practice falls much faster than the lowest order estimate in eq. (78). In a 
recent paper [75], Lappi has shown that the energy density computed in 
this framework, at early times has the form e ~ ln^(l/r); it is finite for any 
r > but is not well defined strictly at r = 0. 

In the discussion up to this point, we have assumed that the color charge 
squared per unit area of the source, g^/J,^, is constant. However, for finite 
nuclei, this is not true and one can define an impact parameter dependent 
Ag, i.e. As{x±). This generalization, in the classical Yang-Mills framework 
described here was discussed previously in Ref. [49] and is given by 

A^g(x^) = A%T^ix^) , (81) 

where Tj^{x±) = J^^dz Py^^g{z,x±) is the nuclear thickness profile, x± is 
the transverse coordinate vector (the reference frame here being the center of 
the nucleus), p^y<,(z, Xpcrp) is the Woods-Saxon nuclear density profile, and 
A^g is the color charge squared per unit area in the center of the nucleus. One 
can use this expression to compute the multiplicity as a function of impact 
parameter in the collision. Then, by using a Glauber model to relate the 
average impact parameter to the average number of participants [76], one 
can obtain the dependence of the multiplicity on the number of participants. 

Previous computations of the centrality dependence of the multiplicity 
and of rapidity distributions were performed in the KLN approach [73,76, 
77]. There however, unlike eq. (81), the saturation scale depends on the 
number of participant nucleons: 

Q^six±) ~ Npart(a;±) , (82) 



with 



B 



(83) 
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In this formula, a^^ is the nucleon-nucleon cross-section, and b the impact 
parameter between the two nuclei. Note that as this form involves the 
thickness functions of both nuclei A and B, it is manifestly not universal - 
in contrast to the definition in eq. 81. For the centrality dependence of the 
multiplicity distributions, the saturation scales defined through eqs. (81) 
or (82) lead to very similar results. This is because the multiplicity, at 
any particular x±, depends on the lesser of the two saturation scales, say 
Qs,A- The dependence on the "non-universal" factor in eq. (83) is then weak 
because, by definition, Tg is large. 

However, the two prescriptions can be distinguished by examining a 
quantity of phenomenological importance, the eccentricity e defined as 

^ J (Px± e{x±) (y2 + x2) ■ 

This quantity is a measure of the asymmetry of the overlap region between 
the two nuclei in collisions at non zero impact parameter. In an ideal hy- 
drodynamical description of heavy ion collisions, a larger initial eccentricity 
may lead to larger elliptic flow [78] than observed, thereby necessitating 
significant viscous effects. Comparisons of model predictions, with different 
initial eccentricities, to data may therefore help constrain the viscosity of 
the Quark Gluon Plasma. In fig. 6, we show results for the eccentricity 
from the kj_ factorized KLN approach with the saturation scale defined as 
in eq. (82) compared to the classical Yang-Mills (CYM) result computed 
with the definition in eq. (81). The KLN Npart definition of Qg leads to the 
largest eccentricity. The universal CYM definition gives smaller values of e 
albeit larger than the traditional parameterization (used in hydrodynamical 
model computations) where the energy density is taken to be proportional 
to the number of participating nucleons. This result is also shown to be 
insensitive to two different choices of the infrared scale m which regulates 
the spatial extent of the Coulomb tails of the gluon distribution. A qual- 
itative explanation of the differences in the eccentricity computed in the 
two approaches is given in Ref. [79] - we refer the reader to the discussion 
there. Also, a further elaboration of the discussion of Refs. [80,81] was very 
recently presented in Ref. [82] -the revised curves our closer to the result in 
Ref. [79]. 

Our discussion thus far of the Glasma has assumed strictly boost invari- 
ant initial conditions on the light cone, of the form specified in eq. (71). 
However, this is clearly an idealization because it requires strict (5-function 
sources as in eq. (62), and that one completely disregards quantum fluctu- 
ations. Because the collision energy is ultra-relativistic and because quan- 
tum fluctuations are suppressed by one power of ag, this was believed to 
be a good approximation. In particular, it was not realized that violations 
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Fig. 6. The eccentricity as a function of impact parameter. The classical field 
CGC result with two different infrared cutoffs m is denoted by GYM [79]. The 
traditional initial eccentricity used in hydrodynamics is a linear combination of 
mostly "Glauber Npart" and a small amount of "Glauber Ncoii". The "KLN" curve 
is the eccentricity obtained from the GGG calculation in Refs. [80,81]. 



of boost invariance lead to a non-Abelian version of the Weibel instabil- 
ity [55] well known in electromagnetic plasmas. To understand the poten- 
tial ramifications of this instability for thermalization, let us first consider 
where the boost invariant results lead us. From eq. (80), it is clear that 
the energy density is far from thermal - in which event, it would decrease 
as r~^/^. The momentum distributions become increasingly anisotropic: 
{k±) Qs and {kz) ~ 1/t. Once the particle-like modes of the classical 
field {k±_ > As) begin to scatter, the occupation number of the field modes 
begins to decrease. How this occurs through scattering was outlined in an 
elegant scenario dubbed "bottom up" by Baier et al [83]. At very early 
times, small elastic scattering of gluons with k±_ ~ Qs dominates and is 
responsible for lowering the gluon occupation number. The Debye mass 
'TT'^ ~ Q1/{Qst) (see fig. 5) sets the scale for these scattering, and the typ- 
ical kz is enhanced by collisions. One obtains kz ~ Qs/ {QstY^^ ■ From this 
dependence, one can estimate that the occupation number of gluons is / < 1 

Q I<2 

for T > as IQs- For proper times greater than these, the classical field 
description becomes less reliable. In the bottom up scenario, soft gluon radi- 
ation from 2 — > 3 scattering processes becomes important at r ~ ^^'^ /Qs- 
The system thermalizes shortly thereafter at rtherm = ois /Qs with a 
temperature T ~ o?J^Qs- The thermalization time scale in this scenario is 
parametrically faster than that obtained by solving the Boltzmann equation 
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for 2^2 processes, which gives Tthcrm exp(as )/Qs [84,85] 

The Debye mass scale is key to the power counting in the bottom up 
scenario. However, as pointed out recently [87] this power counting is af- 
fected by an instability that arises from a change in sign of the Debye 
mass squared for anisotropic momentum distributions [89]. The instability 
is the non-Abelian analog of the Weibel instability [55] in electromagnetic 
plasmas and was discussed previously in the context of QCD plasmas by 
Mrowczynski [90]. One can view the instability, in the configuration space 
of the relevant fields, as the development of specific modes for which the 
effective potential is unbound from below [91]. Detailed simulations in the 
Hard-Loop effective theory in 1 + 1-dimensions [87,91,92] and in 3 + 1- 
dimensions [93,94] have confirmed the existence of this non-Abelian Weibel 
instability. Particle-field simulations of the effects of the instability on ther- 
malization have also been performed recently [95-97]. 

All of these simulations consider the effect of instabilities in systems 
at rest. However, as discussed previously, the Glasma expands into the 
vacuum at nearly the speed of light. Are they seen in the Glasma ? No such 
instabilities were seen in the boost invariant 2 -|- 1-d numerical simulations. 
In the rest of this lecture, we will discuss the consequences of relaxing boost 
invariance in (now) 3 + 1-d numerical simulations of the Glasma fields 
as may be anticipated, non-Abelian Weibel instabilities also arise in the 
Glasma. 

In heavy-ion collisions, the initial conditions on the light cone are never 
exactly boost invariant. Besides the simple kinematic effect of Lorentz con- 
traction at high energies, one also has to take into account quantum fluc- 
tuations at high energies. For instance, as we discussed in the last lecture, 
we will have small quantum fluctiiations at NLO, for each configuration 
of the color sources, which are not boost invariant. Parametrically, from 
the power counting discussed there, quantum fluctuations may be of order 
unity, compared to the leading classical fields which are of order 1/g. 

In the following, we will discuss two simple models of initial conditions 
containing rapidity dependent fluctuations. A more complete theory should 
specify, from first principles, the initial conditions in the boost non-invariant 
case. We will discuss later some recent work in that direction. The only con- 
dition we impose is that these initial conditions satisfy Gauss' law. We con- 
struct these by modifying the boost-invariant initial conditions in eq. (71) 



For recent numerical studies of thermalization due to scattering, see Ref. [86]. 

For a modified "bottom up" scenario, we refer the reader to Ref. [88]. 

This discussion of 3 -|- 1-d numerical simulations is based on work in Refs. [53,54] 
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to 

Ei{0,r],xs_) = 6Ei{r],x±) , 

Er,iO,r],x±) = ig[a{,ai] + 6Er,iri,x±,) , (85) 

while keeping Ai,Arj unchanged. The rapidity dependent perturbations 
5Ei,6Er) are in principle arbitrary, except for the requirement that they 
satisfy the Gauss law. For these initial conditions, it takes the form 

DidEi + D^Er^ = . (86) 

The boost invariance violating perturbations are constructed as follows. 

• We first generate random configurations 5Ei{x±) with 

{5Ei{xi_)5Ej{yj_)) = %5(£c_l - y^} . 



Next, for our first model of rapidity perturbations, we generate a 
Gaussian random function F{ri) with amplitude A 

{F{rj)F{rj')) = A'6{rj-v'). (87) 

For the second model, we also generate a Gaussian random function, 
but subsequently remove high-frequency components of F{r]) 

F{ri) ^ F{rf) = e"^^" e'^''^^ j dr] e^"^ F(r?) , (88) 

where b acts as a "band filter" suppressing the high frequency modes. 
This model is introduced because the white noise Gaussian fluctua- 
tions of the previous model leads to identical amplitudes for all modes. 
As a consequence, the high momentum modes dominate bulk observ- 
ables such as the pressure. The unstable modes we wish to focus on are 
sensitive to infrared modes at early times but their effects are obscured 
by the higher momentum modes from the white noise spectrum. This 
is particularly acute for large violations of boost invariance. There- 
fore, damping these high frequency modes allows us to also study the 
effect of instabilities for larger values of A, or "large seeds" that violate 
boost-invariance. 

For both models, once F{r)) is generated, we obtain for the fluctuation 
fields 

5Ei{r],x±) = d^F{ri)5Ei{x±^) , 

SErjiri, x^) = -F{ri)DiSEi{x^) . (89) 
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These fluctuations, by construction, satisfy Gauss' law. To implement ra- 
pidity fluctuations in the above manner, one requires r > 0. This is a 
consequence of the r, 7] coordinates^'^ and docs not have a physical ori- 
gin. We therefore implement these initial conditions for r = Tinit with 
< Tinit ^ Qs^- Our results are only weakly dependent on the specific 
choice of Tinit- 

The primary gauge invariant obscrvables in simulations of the classical 
Yang-Mills equations are the components of the energy-momentum ten- 
sor [46,49]. We will discuss specifically 

T"^^ + Tyy = 2TT[F^y + E^] , 

Note that the Hamiltonian density is 7^ = tT'^^ = T(r^^ + T^v + t'^T'^). 
These components can be expressed as 

tp^ = I (T^^ + ryy) , tp^ = t^t"^ , (91) 

which correspond to t times the mean transverse and longitudinal pressure, 
respectively. 

When studying the time evolution of rapidity-fluctuations, it is useful 
to introduce Fourier transforms of observables with respect to the rapidity. 
For example, 

P^ (t, u,k^ = d) = j e">'' {P, (t, r?, , (92) 

where denotes averaging over the transverse coordinates {x,y). Apart 
from u = 0, this quantity would be strictly zero in the boost-invariant 
(A = 0) case, while for non-vanishing A and v, P^{y) has a maximum 
amplitude for some specific momentum v. Using a very small but finite vahic 
of A, this maximum amplitude is very much smaller than the corresponding 
amplitude of a typical transverse momentum mode. 

The physical parameters in this study are g^ji (=As; see the discussion 
after cq. (76)), = t:B? , where R is the nuclear radius, A, the initial size 
of the rapidity dependent fluctuations and finally, the band filter h, which as 
discussed previously, we employ only for large values of A. Physical results 
are expressed in terms of the dimensionless combinations fiT and jiL. 
For RHIC collisions of gold nuclei, one has g^nL k, 120; for collisions of lead 

This system of coordinates becomes singular when r ^ 0, as can be seen from the 
fact that the Jacobian for the transformation from Cartesian coordinates vanishes in 
this limit. 
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nuclei at LHC energies, this will be twice larger. The physical properties 
of the spectrum of fluctuations (specified in our simple model here by A 
and b) will presumably be further specified in a complete theory. For our 
present purposes, they will be treated as arbitrary parameters, and results 
presented for a large range in their values. 

Briefly, the lattice parameters in this study, in dimensionless units, are 
(i) A^_|_ and A'^^, the number of lattice sites in the x± and directions 
respectively; (ii) g^^a± and a^, the respective lattice spacings; (iii) TQ/a± 
and 5t, the time at which the simulations are initiated and the stepping 
size respectively. The continuum limit is obtained by holding the physical 
combinations g'^fJ'a±N± = g'^fiL and a^j N^j = fixed, while sending 6t, 
g^Ha± and to zero. For this study, we pick = 1.6 units of rapidity. 
The magnitude of violations of boost invariance, as represented by A, is 
physical and deserves much study. The initial time is chosen to ensure that 
for A = 0, we recover earlier 2 + 1-d results; we set tq = 0.05 a^. Our results 
are insensitive to variations that are a factor of 2 larger or smaller than this 
choice. For further details on the numerical procedure we refer the reader 
to Refs. [53,54]. 

0.0001 

le-05 

le-06 

-i le-07 

"'™ le-08 

1- le-09 

g le-10 
P 

le-11 
le-12 

''^"'^0 500 1000 1500 2000 2500 3000 3500 

2 

g 

Fig. 7. The maximum Fourier mode amplitudes of the mean longitudinal pressure 
Pi = T^T'''' for g'^fiL = 67.9, N± = = 64, N^Ur, = 1.6. Also shown are best fits 
with exp(ar) and exp(a-yr) behavior. The former is clearly ruled out by the data. 

In fig. 7, we plot the maximal value^^ of t'^T^^ = at each time 
step, as a function of g'^fJ-r. The data are for a 64'^ lattice and correspond 
to g^i-iL = 67.9 and = 1.6. The maximal value remains nearly constant 
until ~ 250, beyond which it grows rapidly. A best fit to the functional 
form co+ci exp(c2T'^^) gives C2 = 0.427±0.01 for C3 = 0.5; the coefficients cq. 




What we mean by this is further clarified by fig. 8 and the related discussion there. 
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ci are small numbers proportional to the initial seed. It is clear from fig. 7 
that the form exp(ry^^/zr) is preferred to a fit with an exponential growth 
in r. This eycp{T\/'^^) growth of the unstable soft modes is closely related 
to the mass generated by the highly non-linear dynamics of soft modes in 
the Glasma. As we discussed previously, and showed in fig. 5, a plasmon 
mass Wpi = (jj{k± = 0), is generated. After an initial transient behavior, it 
is of the form 

a;pi = «o5V^p^ • (93) 

with kq = 0.3 ± 0.01 (this parameterization is robust as one approaches the 
continuum limit). The dependence on g^^iL is weak. 

In the finite temperature Hard Thermal Loop (HTL) formalism for 
anisotropic plasmas, the maximal unstable modes of the stress-energy ten- 
sor grow as exp(27r)), where the growth rate 7 satisfies the relation 7 = 
moo/\/2 for maximally anisotropic particle distributions [87]. Here 

where f{p) is the anisotropic single particle distribution of the hard modes. 
It was shown in Ref. [89] that = 3a;pj/2 for both isotropic and aniso- 
tropic plasmas. One therefore obtains 

2 

2^T = \/3 =r = V^hq \/ g^\i t . (95) 

For g^jjiL = 67.9, kq « 0.26 gives the coefficient s/Skq ~ 0.447, which is 
quite close to the value obtained by a fit to the numerical data Ffit = 0.427. 
However, this agreement is misleading because a proper treatment would 
give in the exponent J^cir'7(r') = TthVy/^j with Fth = 2\/3ko- The 
observed growth rate is approximately half of that predicted by directly 
applying the HTL formalism to the Glasma. Despite obvious similarities, 
it is not clear that the equivalence can be expected to hold at this level of 
accuracy. Nevertheless, the similarities in the two frameworks is noteworthy 
as we will discuss now. _ 

In fig. 8 we show the ensemble-averaged rP^ (r, v) for four different sim- 
ulation times. The earliest time {g'^fJiT ~ 22) shows the configuration before 
the instability sets in. At the next time, one sees a bump abovc^ the back- 
ground, corresponding to the distribution of unstable modes. The unstable 
mode with the biggest growth rate (the cusp of the "bumps" in fig. 8) was 
precisely what was used to determine the maximal growth rate F by fitting 
the time dependence of this mode to the form cq + ci exp {T^/g^fir). The 
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Fig. 8. rP^ (t, v) as a function of momentum i^, averaged over 160 initial conditions 
on a 16^ X 2048 lattice with g^^J.L = 22.6 and L,, = 102.4, A ~ 10"". Four different 
simulation times show how the softest modes start growing with an distribution 
reminiscent of results from Hard-Loop calculations [89]. Also indicated are the 
respective values of fmax for three values of g^fir (see text for details). 



two later time snapshots shown in fig. 8 (for g'^fir ~ 156 and g'^fJ-r ~ 288) 
indicate that the growth rate of the unstable modes closely resembles the 
analytic prediction from Hard-Loop calculations [87,89]. 
In fig. 8, 

'^max is the largest mode number that is sensitive to the in- 
stability. Its behavior is shown in fig. 9. From this figure, one observes 




Fig. 9. Left: Time evolution of I'max, on a lattices with g^iiL — 22.7, L,, — 1.6 
and various violations of boost-invariance A. The dashed line represents the linear 
scaling behavior. Right: Time evolution of the maximum amplitude tP^(t, r/). 
When this amplitude reaches a certain size (denoted by the dashed horizontal 
line), J^max starts to grow fast. 
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an underlying trend indicating a linear increase of fmax with approximately 
t'max — 0.06 ^T. For sufficiently small violations of boost-invariance, this 
seems to be fairly independent of the transverse or longitudinal lattice spac- 
ing we have tested. For much larger violations of boost-invariance - or 
sufficiently late times - one observes that i^max deviates strongly from this 
"linear law". In fig. 9 we show that this deviation seems to occur when 
the maximum amplitude of tPl{t, v) reaches a critical size, independent of 
other simulation parameters. This critical value is denoted by a dashed hor- 
izontal line and has the magnitude 3- 10~^ in the dimensionless units plotted 
there. A possible explanation for this behavior is that once the transverse 
magnetic field modes in the Glasma (with small A;^) reach a critical size, the 
corresponding Lorenz force in the longitudinal direction is sufficient to bend 
"particle" (hard gauge mode) trajectories out of the transverse plane into 
the longitudinal direction. This is essentially what happens in electromag- 
netic plasmas. Note however that in electromagnetic plasmas the particle 
modes are the charged fermions.In contrast, the particle modes here are 
the hard ultraviolet transverse modes of the field itself. We will comment 
shortly on how this phenomenon may impact thermalization. 



Fig. 10. Time evolution of the (ensemble-averaged) maximum amplitude of 
TP^(r, I/), for g'^nL = 22.6, = 1.6, N± = 16,32, A = 0.1 a^^ and N^, ranging 
from 32 to 2048. Larger lattices correspond to smaller A. This explains why the 
early-time behavior is not universal for the simulations shown here. 




The saturation seen in the right of fig. 9 is shown clearly in fig. 10 
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where we plot the temporal evolution of the maximum amplitude of the 
ensemble averaged tPl{t,u)^ for lattices with different a^. Early times in 
this figure {g^HT < 200) correspond to the stage when the Weibel instability 
is operative. Interestingly, the simulations show saturation of the growth at 
approximately the same amplitude. These preliminary results are similar 
to the phenomenon of "non-Abelian saturation", found in the context of 
simulations of plasma instabilities in the Hard Loop framework [93,94]. 

In the small seed case, the longitudinal fluctuations carry a tiny fraction 
of the total system energy. In the large seed case, in contrast, for the sim- 
ulations shown here, the initial energy contained in the longitudinal modes 
is ~ 1% of the total system energy. In reality, we expect this fraction to 
be significantly larger. However, this would require us to study the dynam- 
ics on even larger longitudinal lattices than those included in this study to 
ensure that the contributions to the pressure from ultraviolet modes are 
not contaminated by lattice artifacts. In the left fig. 11, we plot as a 
function of r for different lattice spacings o,,. For large (low lattice UV 
cutoff), the longitudinal pressure is consistent with zero; it is clearly finite 
when the lattice UV cutoff is raised. However, the rise saturates as there is 
no notable difference between the simulations for the three smallest values 
of the lattice spacing. At face value, this result suggests that the rise in the 
longitudinal pressure is physical and not a discretization artifact. Clearly, 
further studies on larger transverse lattices are needed to strengthen this 
claim. 

In the right of fig. 11, we investigate the time evolution of the transverse 
pressure and the energy density for (i) a simulation with a low UV cutoff 
(16^ X 32 lattice) and (ii) a simulation with a high UV cutoff (16^ x 2048 
lattice). We observe that the rise in the mean longitudinal pressure ac- 
companies a drop both in the mean transverse pressure and energy density. 
This result is consistent with the previously advocated physical mechanism 
of the Lorenz force bending transverse UV modes (thereby decreasing the 
transverse pressure) into longitudinal UV modes (simultaneously raising the 
longitudinal pressure), thereby pushing the system closer to an isotropic 
state. The energy density depends on the proper time as e ~ ^-1 067^ 
which, while not the free streaming result of e ~ r~^, is also distinct from 
the e ~ r~^/^ required for a locally isotropic system undergoing one di- 
mensional expansion. Furthermore, the time scales (noting that for RHIC 
energies g^n = Ag 2 GeV) are much larger than the time scales of in- 
terest for early thermalization of the Glasma into a QGP. Similar results 



Note here that for the case of "large seed" instabilities, we consider the model with 
the "band filter" b — 0.5, in order to suppress the ultraviolet modes in the initial 
fluctuation. Note further that for larger seeds the instability systematically saturates 
at earlier times, as is clear from the right of fig. 9. 
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Fig. 11. Left: Time evolution of the ensemble and volume averaged longitudinal 
pressure P^, for lattices with g^^L = 22.6, = 1.6, N± = 16,32, A — Cla^^ 
and Njj ranging from 32 to 2048. Note : reduced statistical ensemble of 2 runs for 
A^,; = 4096. Right: Hamiltonian density (_ff (r, ry))^, T(P^(r, 77))^ and r(P^ (r, 77))^, 
for 16^ X 32 and 16^ x 2048 lattices. The energy density is fit to e ~ ^-^-067 
late times. All curves are calculated on lattices with g^^L = 22.6, L,, = 1.6 and 
A = O.lay^ 



were obtained in an analytical model of the late time behavior of expanding 
anisotropic fields in the Hard Loop formalism [98]. 

Nevertheless, these simulations are proof in principle that non-trivial dy- 
namics can take place in the Glasma driving the system towards equilibrium. 
A mode analysis along the lines of that performed recently by Bodeker and 
Rummukainen [99] is required to understand this dynamics in greater detail. 
In particular, it would be useful to understand whether the rapid shift of 
unstable modes to the ultraviolet (as seen in Fig. 9) is due to a turbulent 
Kolmogorov cascade as discussed in Refs. [100,101]. The most important 
task however is to understand from first principles the spectrum of initial 
fluctuations that break boost invariance. A first step in this direction was 
taken in Ref. [102]. This issue is closely related to the NLO computation 
of small fluctuations outlined in Lecture I. At NLO, some of these quantum 
fluctuations are accompanied by large logs in x. Thus for UgY ~ 1, these 
effects are large. To completely understand which contributions from the 
small fluctuations can be absorbed in the evolution of the initial wavefunc- 
tions and to isolate the remainder that contributes to the spectrum of 



Note that the weight functionals describing the distribution of color sources, in this 
case, will have a dependence on rapidity described by the JIMWLK RG equations. 
This implies that even if the underlying classical fields are nearly boost invariant, the 
spectrum of gluons obtained by averaging over all color configurations, are not boost 
invariant. One expects significant deviations from boost invariance when rapidity 
varies by an amount SY ~ 1 /as . 
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initial fluctuations, requires that we demonstrate factorization for inclusive 
multiplicities. This work is in progress [103]. Finally, another interesting 
problem is whether one can match the temporal evolution of Glasma fields 
into kinetic equations at late times. Such a matching was considered pre- 
viously in Ref. [104]. The early time strong field dynamics may however 
modify the power counting assumed in these studies - this possibility is 
also under active investigation [105]. In conclusion, understanding the early 
classical field dynamics of the Glasma and its subsequent thermalization is 
crucial to understand how and when the system thermalizes to form a QGP. 
The phenomenological implications of these studies are significant because 
they influence the initial conditions for hydrodynamic models. One such 
example that we discussed is the initial eccentricity of the QCD matter; its 
magnitude may be relevant for our understanding of just how "perfect" , the 
perfect fluid created at RHIC is. 

Lecture III: Limiting Fragmentation in the CGC framework 

In the flrst two lectures, we discussed the problem of multi-particle pro- 
duction for hadronic collisions where the large x modes are strong sources 
pi.2 ~ 1/5- This is a good model of the dynamics in proton-proton collisions 
at extremely high energies or in heavy ion collisions already at somewhat 
lower energies In eq. (1), one has pi^2/k']_ ~ 1, where k±_ is the typical 
momentum of the partons in the nuclei. As we then discussed in lectures 
I and II, there is no small parameter in the expansion in powers of these 
sources and one has to solve classical equations of motion numerically to 
compute the average inclusive multiplicities for gluon and quark production. 
However, for asymmetric collisions, the most extreme example of which are 
collisions of protons with heavy nuclei, one has a situation where Pp/k\ <^ 1 
and PA/kj_ ~ 1. The other situation where a similar power counting is appli- 
cable is when one probes forward (or backward) rapidities in proton-proton 
or nucleus-nucleus collisions. In these cases, one is probing large x par- 
ton distributions in one hadron (small color charge density - pi/k'j_ <C 1) 
and small x parton distributions in the other (large color charge density - 
P2/k\ ~ 1). In these situations, analytical computations are feasible in the 
CGC framework. In this lecture, we will discuss the phenomenon of limiting 
fragmentation in this framework. 

The hypothesis of limiting fragmentation [106] in high energy hadron- 
hadron collisions states that the pseudo-rapidity distribution ^ (where 
rj' = T] — Ibeam is the pseudo-rapidity shifted by the beam rapidity Ibeam = 
ln(^/s/mp)) becomes independent of the center-of-mass energy i/s in the 



This follows from the large density of color charges in the transverse plane of the 
Lorentz contracted heavy nuclei. 
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region around rj ^ 0, i.e. 

where h is the impact parameter. 

Limiting fragmentation appears to have a wide regime of vahdity. It was 
confirmed experimentally in pp, pA, vryl and AA collisions at high energies 
[107-110]. More recently, the BRAHMS and PHOBOS experiments at the 
Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Labora- 
tory (BNL) performed detailed studies of the pseudo-rapidity distribution 
of the produced charged particles dNch/drj for a wide range (—5.4 < rj < 
5.4) of pseudo-rapidities, and for several center-of-mass energies {y/s^ = 
19.6, 62.4, 130 and 200 GeV) in nucleus-nucleus (Au-Au and Cu-Cu) and 
deuteron-nucleus (d-Au) collisions. Results for pseudo-rapidity distributions 
have also been obtained over a limited kinematic range in pseudo-rapidity 
by the STAR experiment at RHIC [111]. These measurements have opened 
a new and precise window on the limiting fragmentation phenomenon. 

It is worth noting that this scaling is in strong disagreement with boost 
invariant scenarios which predict a fixed fragmentation region and a broad 
central plateau extending with energy. It would therefore be desirable to 
understand the nature of hadronic interactions that lead to limiting frag- 
mentation, and the deviations away from it. In a recent article, Bialas and 
Jezabek [1] , argued that sonic qualitative features of limiting fragmentation 
can be explained in a two-step model involving multiple gluon exchange 
between partons of the colliding hadrons and the subsequent radiation of 
hadronic clusters by the interacting hadrons. Here we will discuss how the 
limiting fragmentation phenomenon arises naturally within the CGC ap- 
proach - we shall address its relation to the Bialas-Jezabek model briefly 
later. 

Inclusive gluon production in proton-nucleus collisions was first com- 
puted in Refs. [113,114], and shown to be k± factorizable in Ref. [115]. In 
Ref. [116], the gluon field produced in pA collisions was computed explic- 
itly in Lorentz gauge df^Af^ = 0. More recently, the gluon field was also 
determined explicitly in the A~^ = light-cone gauge [121]. The inclu- 
sive multiplicity distribution of produced gluons^^ can be expressed in the 
fc^-factorized form as [15,115,116], 

(97) 



Our discussion is based on our paper with A. Stasto [112]. 

For more complicated final states, like quark-antiquark pairs, fcx factorization is 
explicitly violated [117-120]. 
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The formula, as written here, is only valid at zero impact parameter and 
assumes that the nuclei have a uniform density in the transverse plane; the 
functions (p^ g are defined for the entire nucleus. S^^ denotes the transverse 
area of the overlap region between the two nuclei, while 7ri?^ ^ are the total 

transverse area of the nuclei, and Cp = {N^ — l)/2Nc is the Casimir in the 
fundamental representation. 

The longitudinal momentum fractions xi and X2 are defined as 

nip nip 

where ibeam = lii(\/i/mp) is the beam rapidity, nip is the proton mass, and 
Pj^ is the transverse momentum of the produced gluon. The kinematics here 
is the 2 — > 1 eikonal kinematics, which provides the leading contribution to 
gluon production in the CGC picture. 

The functions 0^ and (f)g are obtained from the dipole-nucleus cross- 
sections for nuclei A and B respectively, 

'^-■-^^'^^^ ^ J '^'^^ {U\0)U{X^)))^ , (99) 

where Y = ln(l/x) and where the matrices U are adjoint Wilson lines 
evaluated in the classical color field created by a given partonic configuration 
of the nuclei ^ or in the infinite momentum frame. For a nucleus moving 
in the —z direction, they are defined to be 



Uix±) = T+exp 



+ 00 



(100) 



Here the are the generators of the adjoint representation of SU{Nc) 
and T+ denotes the "time ordering" along the z+ axis. pa{z'^,x±) is a 
certain configuration of the density of color charges in the nucleus under 
consideration, and the expectation value ( • • • ) corresponds to the average 
over these color sources 

As discussed previously, in the McLerran-Venugopalan (MV) model [19- 
21], where no quantum evolution effects are included, the p's have a Gaus- 
sian distribution, with a 2-point correlator given by 

{paiO)Pb{x±)) = p-lSabS^^\x± - y^) , 

where /x^ = A/IttR^ is the color charge squared per unit area. This deter- 
mines (f)j^g completely [115,116], since the 2-point correlator is all we need 
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to know for a Gaussian distribution. We will shortly discuss the small x 
quantum evolution of the correlator on the r.h.s. of eq. (99). 

These distributions 4>^ g , albeit very similar to the canonical unintc- 
grated gluon distributions in the hadrons, should not be confused with the 
latter [115,116]. However, at large k± {k± ^ Qs), they coincide with the 
usual unintegrated gluon distribution. Note that the unintegrated gluon dis- 
tribution here is defined such that the proton gluon distribution, to leading 
order satisfies 



1 f'^ 

xGp{x,Q'^) = J dl']_(()p{x,l±) 



From eq. (97), it is easy to see how limiting fragmentation emerges in 
the limit where X2 <C xi. In this situation, the typical transverse mo- 
mentum k± in the projectile at large xi is much smaller than the typical 
transverse momentum — k± [ in the other projectile, because these are 
controlled by saturation scales evaluated respectively at xi and at X2 respec- 
tively. Therefore, at sufficiently high energies, it is legitimate to approxi- 
mate <Pb{x2, \p± — k±\) by <pg{x2,p±)- Integrating the gluon distribution 
over pj^, we obtain 

^ = 2.^CA^Rl)i.Rl) J -^<i>six^^P-) J ^2 

= ^<^.(^i,fc±)^ ^^15(^1, m'^q^m). (101) 



This expression is nearly independent of X2 and therefore depends only 
weakly on on y — ybeam- To obtain the second line in the above expression, 
we have used eq. (99) and the fact that the Wilson line C/ is a unitary 
matrix. Therefore, details of the evolution are unimportant for limiting 
fragmentation, only the requirement that the evolution equation preserves 
unitarity. The residual dependence on X2 comes from the the upper limit 
~ Qf{x2) of the integral in the second line. This ensures the applicability of 
the approximation that led to the expression in the second line above. The 
integral over k± gives the integrated gluon distribution in the projectile, 
evaluated at a resolution scale of the order of the saturation scale of the 
target. Therefore, the residual dependence on y -|- Ibeam arises only via the 
scale dependence of the gluon distribution of the projectile. This residual 
dependence on y + Ibcam is very weak at large xi because it is the regime 
where Bjorkcn scaling is observed. 

The formula in eq. (101) was used previously in Ref. [122]. The nu- 
clear gluon distribution here is determined by global fits to deeply inelastic 
scattering and Drell-Yan data. We note that the glue at large x is very 
poorly constrained at present [123]. The approach of Bialas and Jezabek [1] 
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also amounts to using a similar formula, although convoluted with a frag- 
mentation function (see eqs. (1), (4) and (5) of [1] - in addition, both the 
parton distribution and the fragmentation function are assumed to be scale 
independent in this approach). We will discuss the effect of fragmentation 
functions later in our discussion. 

Though limiting fragmentation can be simply understood as a conse- 
quence of unitarity in the high energy limit, what may be more compelling 
are observed deviations from limiting fragmentation and how these vary 
with energy. We will now see whether deviations from limiting fragmenta- 
tion can be understood from the renormalization group (RG) evolution of 
the unintegrated gluon distributions in eq. 97. In particular, we study the 
RG evolution of these distributions given by the Balitsky-Kovchegov (BK) 
equation [33,34]. 

The BK equation is a non-linear evolution equation^*^ in rapidity Y = 
ln(l/a;) for the forward scattering amplitude T{r±, Y) of a quark-antiquark 
dipole of size r± scattering off a target in the limit of very high center-of- 
mass energy -s/s where T is defined as: 

T(rx, Y) = l- i-TV (^U^{0)U{r^))^ . (102) 

Here U is the corresponding Wilson line for the scattering of a quark- 
anti-quark dipole in the fundamental representation. The correlators U 
in eq. (99), which are in the adjoint representation are Wilson lines for the 
scattering of a gluon dipole on the same target instead. The BK equation 
captures essential features of high energy scattering. When r± <^ ^/Qs, 
one has color transparency; for r± 3> ^/Qsi the amplitude T — > 1, and one 
obtains gluon saturation It is therefore an excellent model to study both 
limiting fragmentation as well as deviations from it. 

It is convenient to express the BK equation in momentum space in terms 
of the Bessel-Fourier transform of the amplitude 

+00 

f{k^,Y)= [ —Mk±r^)Tir^,Y) . (103) 


One obtains 

=asiK^f){k^,Y)-asf\k^,Y) , (104) 



It is equivalent to the corresponding JIMWLK equation [25-32] of the Color Glass 
Condensate, in a mean field (large Nc and large A) approximation where higher order 
dipole correlators are neglected. 

Qs is defined in terms of the requirement that T = 1/2 for r± = 2/Qs. 
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where we denote = agNc/Tr. The operator K is the weU known BFKL 
kernel in momentum space [10,11]. 

In the large Nc and large A limit, the correlators of Wilson lines in the 
fundamental and adjoint representations are simply related: 

i- Tr([/(0)C/t(r^))^ = N,[1- T{r^, Y)f . 

One can therefore express the unintegrated gluon distribution in eq. (99) in 
terms of T as 

</'A,s(2;,fc±) = J r±dr± Jo{k_Lr_L) [1 - ^ (r_L, ln(l/x))J . 

(105) 

In Ref. [112] we solved the BK equation numerically, in both fixed and 
running coupling cases, in order to investigate limiting fragmentation in 
hadronic collisions The results are shown in fig. 12. The solid line is 




Fig. 12. The unintegrated gluon distribution from (i) correlators in the adjoint 
representation (eq.105) (dashed lines) and (ii) correlators in the fundamental rep- 
resentation - see text (solid lines). 

the result obtained for the unintegrated distribution corresponding to cor- 
relators in the fundamental representation, i.e. proportional to the Fourier 
transform of 1 — T instead of that of Nc{l — T)^ in eq. (105). 

The BK equation was solved numerically previously for various studies [124-128] 
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Our results for limiting fragmentation are obtained through the following 
procedure: 

• One first solves the BK equation in eq. (104) to obtain (via eq. (103)) 
eq. (105) for the unintegrated distributions 4>{x,k±). The solution is 
performed for x < xq = 10~^ with the initial condition (f>{xo,k±), 
given by the McLerran-Venugopalan model [19-21] with a fixed initial 
value of the saturation scale Qfixo). For a gold nucleus, extrapola- 
tions from HERA and estimates from fits to RHIC data suggest that 
{Qf{x()))'^ ~ 2GeV^. The saturation scale in the proton is taken to 
be Q^{xo) = ((^^(xo))^ (197/^)^/^. For comparison, we also consid- 
ered initial conditions from the Golec-Biernat and Wusthoff (GBW) 
model [129]. The values of Qf were varied in this study to obtain best 
fits to the data. 

• We used the ansatz 

'^^^'^^) = (r^) (?) 

in order to extrapolate our results to larger values of x > xq, where 
the parameter /3 = 4 is fixed by QCD counting rules. 

• The resulting expressions are substituted in eq. (97) to determine ra- 
pidity distribution of the produced gluons. The pseudo-rapidity dis- 
tributions are determined by multiplying eq. (97) with the Jacobian 
for the transformation from y to rj. This transformation requires one 

to specify an infrared mass, which is also the mass chosen to regulate 
the (logarithmic) infrared sensitivity of the rapidity distributions. 

For further details on how the results are obtained, we refer the reader to 
Ref. [112]. 

In figure 13, we plot the pseudo-rapidity distributions of the charged 

particles produced in nuclcon-nuclcon collisions for center of mass energies 
ranging from 53GeV to 900 GeV. The computations were performed for 
input distributions (for BK evolution) at xq = 0.01 from the the GBW and 
MV models. The normalization is a free parameter which is fitted at one 
energy. Plots on the left of figure 13 are obtained for Aq = 0.15 (the free 
parameter in the large x extrapolation) whereas the right plots are for Aq = 
0.0. The different values of A,c = 4.88as arc obtained for different values of 
as as inputs to the BK equation. While these values of ag might appear 
small, they can be motivated as follows. The amplitude has the growth 
rate Abk = 2.77 As/4.88 0.57 Ag. Thus As = 0.46, which gives reasonable 
fits (more on this in the next paragraph) to the pp data for the MV initial 
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Fig. 13. Pseudorapidity distributions dN/drj for charged particles from nucleon- 
nucleon collisions at UA5 energies [107] Vs = 53, 200, 546, 900 GeV and PHOBOS 
data [110] at ^/s — 200 GeV. Upper plots: initial distribution from the MV model, 
lower plots: initial distribution from the GBW model. Left panels: Aq = 0.15, 
right panels Aq = 0.0. 



conditions, corresponds to Abk = 0.28. Thus a small value of in fixed 
coupling computations "mimics" the value for the energy dependence of the 
amplitude in next-to-leading order resummed BK computations [39] and in 
empirical dipole model comparisons to the HERA data [129]. 

Our computations are extremely sensitive to the extrapolation prescrip- 
tion to large x. This is not a surprise as the wave-function of the projectile 
is probed at fairly large values of xi. From our analysis, we see that the 
data naively favors a non-zero value for Aq- The value Aq = results in 
distributions which, in both the MV and GBW cases, give reasonable fits 
(albeit with different normalizations) at lower energies but systematically 
become harder relative to the data as the energy is increased. To fit the 
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data in the MV model up to the highest UA5 energies, a lower value of As 
than that in the GBW model is required. This is related to the fact that 
MV model has tails which extend to larger values in k± than in the GBW 
model. As the energy is increased, the typical k± ~ Qs{x2) does as well. 
We will return to this point shortly. 




Fig. 14. The same as figure 13 but plotted versus rj' ~ rj — Tbcam to illustrate the 
region of limiting fragmentation. 

In figure 14 the same distributions are shown as a function of the rj' = 
V — Yheam- The Calculations for Aq = 0.15 are consistent with scaling in 
the limiting fragmentation region. There is a slight discrepancy between 
the calculations and the data in the mid-rapidity region. This discrepancy 
may be a hint that one is seeing violations of k± factorization in this regime 
because k± factorization becomes less reliable the further one is from the 
dilute-dense kinematics of the fragmentation regions [45,130]. This discrep- 
ancy should grow with increasing energy. However, our parameters are not 
sufficiently constrained that a conclusive statement can be made. For in- 
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stance, as we mentioned previously, there is a sensitivity to the infrared 
mass chosen in the Jacobian of the transformation from y to rj. This is 
discussed further in Ref. [112]. 

In figure 15 we show the extrapolation to higher energies, in particular 
the LHC range of energies for the calculation with the GBW input. We 
observed previously that the MV initial distribution, when evolved to these 
higher energies, gives a rapidity distribution which is very flat in the range 
—5 < 7/ < 5. We noted that this is because the average transverse momen- 
tum grows with the energy giving a significant contribution from the high k± 
tail of the distribution in the MV input at .xq. The effect of fragmentation 
functions on softening the spectra in the limiting fragmentation region can 
be simply understood by the following qualitative argument. The inclusive 
hadron distribution can be expressed as 



Dg_^Jz = ^,^i''] , (106) 



d?p^dy J^^.^^ z d'^q^dy ^'^^ \ q_L 

where Dg_^h{z, fi^) is the fragmentation function denoting the probability, 
at the scale that a gluon fragments into a hadron carrying a fraction 
z of its transverse momentum. For simplicity, we only consider here the 
probability for gluons fragmenting into the hadron. The lower limit of the 
integral can be determined from the kinematic requirement that xi,2 < 1 - 
we obtain, 

2min = — e^-^b— . (107) 

If 

■^min were zero, the effect of including fragmentation effects would simply 
be to multiply eq. (106) by an overall constant factor. At lower energies, 
the typical value of q± is small for a fixed y — ibeam; the value of z^am is 
quite low. However, as the center of mass energy is increased, the typical 
g_L value grows slowly with the energy. This has the effect of raising Zmm for 
a fixed y — Ibeam; thereby lowering the value of the multiplicity in eq.(106) 
for that y — Ibcam- Note further that eq. (107) suggests that there is a 
kinematic bound on q_\_ as a function of y — Ibeam ~ only very soft gluons 
can contribute to the inclusive multiplicity. 

In figure 16 we display the p_L distributions obtained from the MV input 
compare to the UAl data [131]. We compare the calculation with and 
without the fragmentation function. The fragmentation function has been 
taken from [132]. Clearly the "bare" MV model does not describe the data at 
large fc_L because it does not include fragmentation function effects which, as 
discussed, make the spectrum steeper. In contrast, because the k± spectrum 
of the GBW model dies exponentially at large this "unphysical" k±_ 
behavior mimics the effect of fragmentation functions - see figure 16. Hence 
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Fig. 15. Predictions for higher center-of-mass energies for proton-proton coUisions: 
\/s — 2, 8, 14TeV for GBW input modeL The parameter in the large x extrapola- 
tion was set to Aq = 0.15. 



extrapolations of this model, as shown in figure 18 give a more reasonable 
looking result. Similar conclusions were reached previously in Ref. [133]. 

We next compute the pseudo-rapidity distribution in deuteron-gold col- 
lisions. In figure 17 we show the result for the calculation compared with the 
dA data [135]. The unintegrated gluons were extracted from the pp and AA 
data. The overall shape of the distribution matches well on the deuteron 
side with the minimum-bias data. The disagreement on the nuclear frag- 
mentation side is easy to understand since, as mentioned earlier, it requires 
a better implementation of nuclear geometry effects. Similar conclusions 
were reached in Ref. [134] in their comparisons to the RHIC deuteron-gold 
data. 

We now turn to to nucleus-nucleus collisions. In figure 18 we present fits 
to data on the pseudo-rapidity distributions in gold-gold collisions from the 
PHOBOS, BRAHMS and STAR collaborations. The data [110] are for ^ = 
19.6, 130, 200 GeV and the BRAHMS data [109] are for ^ = 130, 200 GeV. 
A reasonable description of limiting fragmentation is achieved in this case 
as well. One again has discrepancies in the central rapidity region as in the 
pp case. We find that values of Q'^ « 1.3 GeV for the saturation scale give 
the best fits. This value is consistent with the other estimates discussed 
previously [48,72,77]. Apparently the gold-gold data are better described 
by the calculations which have Aq = 0.0. This might be related to the 
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Fig. 16. p± distribution from eq. (97) with MV (full squares) and GBW (full 
triangles) initial conditions. The MV initial condition - with the fragmentation 
function included ~ is denoted by the open squares. The distribution is averaged 
over the rapidity region y = 0.0 — 2.5, to compare with data (in 200 GeV/nucleon 
proton-antiproton collisions in the same pseudo-rapidity range) on charged hadron 
p± distributions from the UAl collaboration: full circles. 

difference in the large x distributions in the proton and nucleus. Further, 
slightly higher values of are preferred to the pp case. This variation of 
parameters from AA to pp case might be also connected with the fact that 
in our approach the impact parameter is integrated out thereby averaging 
over details of the nuclear geometry. 

In fig. 19 we show the extrapolation of two calculations to higher energy 
^/s = 5500 GeV. We note that the calculation within the MV model gives 
results which would violate the scaling in the limiting fragmentation region 
by approximately 20% at larger y — ibeam- This violation is due partly to the 
effect of fragmentation functions discussed previously and partly to the fact 
that the integrated parton distributions from the MV model do not obey 
Bjorken scaling at large values of x. In the latter case, the violations are 
proportional to ln{Q'^{x2)) as discussed previously. The effects of the former 
are simulated by the GBW model - the extrapolation of which, to higher 
energies, is shown by the dashed line. The band separating the two there- 
fore suggests the systematic uncertainity in such an extrapolation coming 
from (i) the choice of initial conditions and (ii) the effects of fragmentation 
functions which are also uncertain at lower transverse momenta. 
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Fig. 17. Comparison of theory (with MV input) to minimum bias deuteron-gold 
data at RHIC [135]. 

To summarize the discussion in this lecture, we studied the phenomenon 
of hmiting fragmentation in the Color Glass Condensate framework. In 
the dilute-dense (projectile-target) kinematics of the fragmentation regions, 
one can derive (in this framework) an expression for inclusive gluon distri- 
butions which is k± factorizable into the product of "unintegrated" gluon 
distributions in the projectile and target. From the general formula for 
gluon production (eq. (97)), limiting fragmentation is a consequence of two 
factors: 

• Unitarity of the U matrices which appear in the definition of the un- 
integrated gluon distribution in eq. (99). 

• Bjorken scaling at large xi, namely, the fact that the integrated gluon 
distribution at large x, depends only on xi and not on the scale Qsix2)- 
(The residual scale dependence consequently leads to the dependence 
on the total center-of-mass energy.) 

Deviations from the limiting fragmentation curve at experimentally accessi- 
ble energies are very interesting because they can potentially teach us about 
how parton distributions evolve at high energies. In the CGC framework, 
the Balitsky-Kovchegov equation determines the evolution of the uninte- 
grated parton distributions with energy from an initial scale in x chosen 
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Fig. 18. Pseudorapidity distributions normalized by the number of participants for 
charged hadrons in gold-gold collisions from the PHOBOS collaboration at energies 
19.6, 130, 200 GeV (filled triangles , squares and circles), BRAHMS collaboration 
at energies 130, 200 GeV (open squares and circles) . The data from the STAR 
collaboration at energy 62.4 GeV (open triangles) are not visible on this plot but 
can be seen more clearly in fig. 19. Upper solid line: initial distributions from the 
MV model; lower solid line: initial distributions from the GBW model. 



here to be xq = 0.01. This choice of scale is inspired by model comparisons 
to the HERA data. 

We compared our results to data on limiting fragmentation from pp 
collisions at various experimental facilities over a wide range of collider 
energies, and to collider data from RHIC for deuteron-gold and gold-gold 
collisions. We obtained results for two different models of initial conditions 
at X > Xq; the McLerran-Venugopalan model (MV) and the Golec-Biernat- 
Wusthoff (GBW) model. 

We found reasonable agreement for this wide range of collider data for 
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Fig. 19. Extrapolation of calculations for gold-gold collisions shown in fig. 18 to the 
LHC energy ^/s = 5500 GeV/ nuclcon. For comparison, the same data at lower 
energies are shown. (See fig. 18.) Dashed line - GBW input Ao = 0, As = 0.46, 
sohd line - MV input with Ao 0, A^ 0.46. 

the limited set of parameters and made predictions which can be tested 
in proton-proton and nucleus-nucleus collisions at the LHC. Clearly these 
results can be fine tuned by introducing further details about nuclear ge- 
ometry. More parameters are introduced, however there is more data for 
different centrality cuts - we leave these detailed comparisons for future 
studies. In addition, an important effect, which improves agreement with 
data, is to account for the fragmentation of gluons in hadrons. In particular, 
the MV model, which has the right leading order large k±_ behavior at the 
partonic level, but no fragmentation effects, is much harder than the data. 
The latter falls as a much higher power of k^. As rapidity distributions at 
higher energies are more sensitive to larger k± , we expect this discrepancy to 
show up in our studies of limiting fragmentation and indeed it does. Taking 
this into account leads to more plausible extrapolations of fits of existing 
data to LHC energies. 
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